ICA-based spatiotemporal approach for single-trial analysis of
postmovement MEG beta synchronization
夞
Po-Lei Lee,
aYu-Te Wu,
a,bLi-Fen Chen,
a,cYong-Sheng Chen,
a,dChou-Ming Cheng,
aTzu-Chen Yeh,
a,eLow-Tone Ho,
a,eMau-Song Chang,
a,eand Jen-Chuen Hsieh
a,e,f,g,*
aLaboratory of Integrated Brain Research, Department of Medical Research and Education, Taipei Veterans General Hospital, Taipei, Taiwan bInstitute of Radiological Sciences, National Yang-Ming University, Taipei, Taiwan
cCenter for Neuroscience, National Yang-Ming University, Taipei, Taiwan
dDepartment of Computer Science and Information Engineering, National Chaio Tung University, Hsinchu, Taiwan eFaculty of Medicine, School of Medicine, National Yang-Ming University, Taipei, Taiwan
fInstitute of Neuroscience, School of Life Science, National Yang-Ming University, Taipei, Taiwan
gInstitute of Health Informatics and Decision Making, School of Medicine, National Yang-Ming University, Taipei, Taiwan Received 8 April 2003; revised 3 July 2003; accepted 14 July 2003
Abstract
The extraction of event-related oscillatory neuromagnetic activities from single-trial measurement is challenging due to the non-phase-locked nature and variability from trial to trial. The present study presents a method based on independent component analysis (ICA) and the use of a template-based correlation approach to extract Rolandic beta rhythm from magnetoencephalographic (MEG) measurements of right finger lifting. A single trial recording was decomposed into a set of coupled temporal independent components and corresponding spatial maps using ICA and the reactive beta frequency band for each trial identified using a two-spectrum comparison between the postmovement interval and a reference period. Task-related components survived dual criteria of high correlation with both the temporal and the spatial templates with an acceptance rate of about 80%. Phase and amplitude information for noise-free MEG beta activities were preserved not only for optimal calculation of beta rebound (event-related synchronization) but also for profound penetration into subtle dynamics across trials. Given the high signal-to-noise ratio (SNR) of this method, various methods of source estimation were used on reconstructed single-trial data and the source loci coherently anchored in the vicinity of the primary motor area. This method promises the possibility of a window into the intricate brain dynamics of motor control mechanisms and the cortical pathophysiology of movement disorder on a trial-by-trial basis.
© 2003 Elsevier Inc. All rights reserved.
Keywords: Rolandic rhythm; Motor cortex; Single-trial; Magnetoencephalography; Event-related synchronization; Independent component analysis (ICA)
Introduction
The human brain is a dynamic system that frequently changes functional mode (Lopes da Silva, 1991, 1996). Spatiotemporal analysis of brain activities with regard to
distinct spatial locations and frequency bands reveals task-specific brain activation which changes in a fraction of a second (Jensen and Vanni, 2002). At rest, Rolandic EEG and MEG rhythms are dominated by rhythmic activity around 10 (alpha band) and 20 (beta band) Hz. Electrocor-ticographic (Pfurtscheller et al., 1994) and neuromagnetic
recordings have shown that the ⬃20-Hz rhythm mainly
originates in the anterior bank of the central sulcus while the
⬃10-Hz rhythm is concentrated predominantly in the
post-central cortex (Pfurtscheller and Lopes da Silva, 1999). These two frequency components appear to have different
functional roles, with the ⬃20-Hz rhythm being more
夞The authors are affiliated with Brain Research Center, University
System of Taiwan.
* Corresponding author. Laboratory of Integrated Brain Research, De-partment of Medical Research and Education, Taipei Veterans General Hospital, No.201, Sect.2, Shih-Pai Rd., Taipei 112 Taiwan. Fax: ⫹886-2-28745182
E-mail address: [email protected] (J.-C. Hsieh).
NeuroImage 20 (2003) 2010 –2030 www.elsevier.com/locate/ynimg
1053-8119/$ – see front matter © 2003 Elsevier Inc. All rights reserved. doi:10.1016/j.neuroimage.2003.07.024
closely connected to movements and their termination and the ⬃10-Hz component behaving more like a classical “idling” rhythm (Salmelin et al., 1995). Voluntary move-ment is composed of three phases: planning, execution, and recovery (Pfurtscheller et al., 1998a). It has been suggested that localized event-related alpha desynchronization (ERD) upon movement can be viewed as an EEG/MEG correlate of an activated cortical sensorimotor network, servicing plan-ning, and execution, while beta event-related synchroniza-tion (ERS) may reflect deactivasynchroniza-tion/inhibisynchroniza-tion during the
recovery phase in the underlying cortical network
(Pfurtscheller et al., 1996).
Movement-related ERD and ERS have been used as probes to study neurophysiology in normal brains and pathophysiology in the diseased (Tamas et al., 2003). It has been reported that the diagnostic features of patients with Parkinson’s disease, in comparison with controls, are a slowing and suppression of the postmovement beta ERS independent of the amount of beta activity in the reference period (Pfurtscheller et al., 1998a). These findings imply that slowed and reduced recovery after the motor act
im-pedes cortical preparation of the next movement
(Pfurtscheller et al., 1996). Patients with Unverricht–Lund-borg-type myoclonic epilepsy demonstrate little rebound of
beta activities contingent upon median nerve stimulation (Silen et al., 2000). The diminished beta ERS indicates that the myoclonic patients have sustained motor cortex reactiv-ity which can be attributed to impaired cortical inhibition (Pfurtscheller and Lopes da Silva, 1999).
ERD and ERS activities are time-locked, but not
phase-locked, to external stimuli or tasks (Andrew and
Pfurtscheller, 1995; Kalcher and Pfurtscheller, 1995; Pfurtscheller and Lopes da Silva, 1999). Existing methods for extraction of ERD/ERS signals essentially measure power or amplitude changes of corresponding frequency bands as derived from the average of dozens or hundreds or trials. The band power method squares and averages filtered brain signals within a selected frequency band (Pfurtscheller and Aranibar, 1977), and an intertrial variance method to remove the phase-locked portion in the band power method was reported by Klimesch et al. (1998). Likewise, autore-gressive and spectral decomposition methods have been used to extract significant frequency components in rhyth-mic signals (Florian and Pfurtscheller, 1995). Salmelin’s temporal–spectral evolution method rectifies and averages filtered MEG signals (Salmelin et al., 1995). To increase the temporal resolution of the ERD/ERS technique, Clochon et al. (1996) proposed an amplitude modulation (AM) method based on the Hilbert transform to detect the envelope of filtered signals by squaring and summing their real and imaginary parts. All these approaches presume stereotypical frequency and temporal characteristics across trials and re-quire an average of many trials for the ERD/ERS using a preset frequency filter and time window to preprocess every trial. However, non-phase-locked rhythmic signals can vary from trial to trial contingent upon variations in a subject’s performance and state, which may be linked to fluctuations in expectation, attention, arousal, and task strategy (Bas-tiaansen et al., 1999, 2001; Earle, 1988; Haig et al., 1995; Hoffman et al., 1991; Yabe et al., 1993). Since trial-to-trial variability in amplitude, latencies, or scalp distribution might carry important information on cognitive and physi-ological states (Jung et al., 2001), a method that permits the extraction and analysis of the oscillatory signal on a single-trial base is crucial for the study of subtle brain dynamics. Furthermore, such a method should require fewer trials for analysis and hence shorter experiment time, which is ben-eficial for patients with impairment of motor and/or cogni-tive performance (Muller-Gerking et al., 1999).
Single-trial multichannel EEG analysis has been devel-oped for time-locked, phase-locked, evoked brain activities (Jung et al., 2001; Tang et al., 2002). However, approaches to single-trial movement-related oscillatory changes are less explored. ICA, a data-driven method for multivariate data analysis, has been used to reveal temporally independent neuronal activities of EEG measurements (Jung et al., 2001; Makeig et al, 1997; McKeown et al., 1998), MEG measure-ments (Wu et al., 2002, 2003; Tang et al., 2002), fMRI (Duann et al., 2002; McKeown et al., 1998), and recently perfusion MRI (Kao et al., 2003). The present study
pro-Fig. 1. Determination of task-specific frequency band using two 1-s am-plitude spectra. (a) R, represents the reference period from⫺4 to ⫺3 s preceding onset of movement; P, postmovement duration from 0.8 to 1.8 s after onset of movement. (b) Two spectra computed over the reference (R) and postmovement periods (P), respectively. (c) The task-specific fre-quency band for beta-band VAMW is defined as the one where the difference between two spectra exceeds the 95% confidence level.
poses a new approach using ICA and the Hilbert transfor-mation for the single-trial detection of movement-related beta rhythmic activity during a self-paced right finger lifting task. This study focuses on beta activity and beta ERS, centered around 20 Hz, because it has been demonstrated that the movement-related short bursts of beta oscillation have higher task and movement specificity than alpha ERD (Pfurtscheller and Aranibar, 1979; Pfurtscheller et al., 1996).
Since brain oscillation may be expressed alone in a specific frequency band independent of artifacts (Ermer et al., 2000; Lins et al., 1993a, 1993b; Mosher et al., 1992), ICA is applied to transform brain signals across all channels (in a single trial) into mutually independent components by means of an unmixing matrix in which each column repre-sents a spatial map tailoring the weights of the correspond-ing temporal component at each MEG sensor. The spatial maps and temporal waveforms of decomposed independent components are categorized into task-related and task-un-related groups respectively, based on temporal and spatial characteristics. This temporal template is the grand average
of hundreds of vector-norm envelopes of the bandpass-filtered, single-trial MEG measurements obtained from right index finger lifting. The spatial template can be derived from the spatial distribution at beta rebound activity either from the grand average of the generation group (for signal extraction) or from each individual (for verification). Cor-relations between the temporal template and component waveforms, as well as between the spatial template and spatial maps, are computed, and coupled component wave-forms and spatial maps that conjointly survive with high correlation values are taken as task-related information and subjected to data reconstruction. In this way the phase and amplitude information of noise-free MEG beta activities can be preserved for profound studies of temporal and spectral variation across trials. Due to the high SNR in beta activities extracted through ICA, trial-specific reactive frequency ranges can be determined by means of the comparisons of two short time spectra between the reference and postmove-ment periods. Beta reactivity per single trial can be quanti-fied using the amplitude modulation (AM) method (Clochon et al., 1996), and significant epochs can be determined using
Fig. 2. Creation of common temporal and spatial templates. (a) The common temporal template, VAMWtemplate, is created by averaging VAMWs (500 trials, 100 trials for each subject, 5 subjects pooled). Event-related beta modulation is defined as the amplitude difference between the mean amplitude of baseline activity (⫺2.5 to ⫺2 s) and maximum amplitude in the postmovement interval (0.8 to 1.8 s). (b) The common spatial template is the average of the topographical distributions of event-related beta modulations of five subjects from model generation group. Only the half the spatial map (unshaded) contralateral to the side of finger lifting is used as the spatial template.
a nonparametric sign test (Brovelli et al., 2002). Source estimation and localization techniques can be successfully applied to single-trial epoch to estimate the source locations of beta modulation.
The current study presents: (1) a novel ICA-based spa-tiotemporal approach for single-trial analysis of event-re-lated beta oscillatory modulations with a high extraction rate; (2) the prospect of trial-specific frequency bandpass filtering that takes into account subtle trial-by-trial brain dynamics; (3) the feasibility of using sophisticated source estimation/localization methods demanding high SNR on single trial data; and (4) a common template approach permitting an effective alternative in cases where lengthy
procedures cannot be endured by participants or in clinical settings where patients have attention problems or are inca-pable of sustaining long experiments.
Materials and methods
Subjects and task
The present study examined six healthy right-handed subjects (gender balanced), ages 24 –30 years. Five subjects were used in the model generation group, and MEG data from the last subject were used for validation. Subjects
Fig. 3. Examples of IC selection and signal reconstruction procedure. (a) Spatial maps, IC waveforms, Fourier spectra of IC waveforms, and VAMWICs of five ICs obtained from one single epoch by ICA. Only ICs fulfilling the dual criteria are selected for signal reconstruction. For example, IC 3 meets the dual criteria (underscored in red): (i) correlation value between spatial map and spatial template is 0.84 (rank⫽ 97%, Z ⫽ 1.89, P ⫽ 0.03); (ii) correlation value between 16 –20 Hz VAMWICand 20 –24 Hz VAMWICwith VAMWtemplateis 0.8 (rank⫽ 99%, Z ⫽ 3.08, P ⫽ 0.01) and 0.78 (rank ⫽ 97.8%, Z ⫽ 2.85, P⫽ 0.022), respectively. (b) Noise identification and removal. The deselected IC 2 in Fig. 2a may emanate from background noise since it resembles the
IC 1 extracted from empty room measurement. (c) The impact of including task-unrelated IC into signal reconstruction. This figure illustrates (a different trial from Fig. 3a) that inclusion of task-unrelated IC (IC 9) with a high spatial correlation (correlation value⫽ 0.61, rank ⫽ 95.2%, Z ⫽ 1.67, P ⫽ 0.048) but poor temporal correlation (correlation value⫽ 0.28, rank ⫽ 13%, Z ⫽ 0.34, P ⫽ 0.87) causes deterioration in the beta BR from 28.9 fT/cm (arrows and trace in red; IC 9 eliminated from reconstruction) to 18.6 fT/cm (arrows and trace in blue; IC9 included for reconstruction).
Fig. 4. Example of single-trial epoch selection based on a nonparametric sign test. Single-trial VAMWrecon_maxs of reconstructed data are examined through a nonparametric sign test. ZIOI(i) is the Z value of the ith trial, NIOI⫹(i) is the number of data points in postmovement IOI with values larger than the median of baseline activities of the ith trial, and NIOIis the total number of time points in postmovement IOI. Only epochs showing significant increase of beta activities are chosen for further analysis. The first trial with a ZIOIscore equal to⫺4.53 is marked as an insignificant epoch and eliminated from further analysis.
performed self-paced lifting of the right index finger ap-proximately once every 8 s. Subjects were trained to per-form the movement briskly for a duration of 200 to 300 ms, as monitored by surface electromyogram (EMG) on exten-sor digitorum communis, with a range of finger movement around 35– 40°, while keeping their eyes open in order to suppress the occipital alpha rhythm. In addition, somato-sensory evoked fields (SEFs) for right median nerve
stim-ulation were measured to locate the primary sensorimotor area (SMI) in each subject as part of the procedure for the generation of a temporal template (see below).
Data recording
Cortical magnetic signals were recorded with a 306-channel (102 sensor unit) whole-head neuromagnetometer
Fig. 5 (continued)
(bandpass, 0.05–250 Hz; digitized at 1 kHz; Vectorview; Neuromag Ltd., Helsinki, Finland) with subjects in sitting position. Each sensor unit was composed of a pair of planar gradiometers and a magnetometer. The magnetometer mea-sured magnetic flux (Bz), normal to the sensor unit, while
the gradiometers measured two tangential derivatives of Bz
(⭸Bz/⭸x and ⭸Bz/⭸y, mutually orthogonal). Only magnetic
signals measured by the gradiometers were used in this study. Bipolar horizontal and vertical electrooculograms (EOG) were recorded using electrodes placed below and above the left eye and at the bilateral outer canthi to monitor eye movement and blinks. The exact position of the head
Fig. 6. Sensor-array display of VAMWrecons and VAMWs. (a) One example of ICA single-trial VAMWrecons of all sensor sites in subject I. The single-trial result shows only left sensorimotor area dominance of event-related activities, as the present study focuses on the area contralateral to movement side and only the left spatial template is used. The dashed trapezoid marks the nine SMI vicinal sensor sites and the VAMWrecon_maxis marked with the red circle. (b) Trial-specific frequency band used for VAMWrecons calculation in Fig. 6a. (c) VAMWs obtained from the conventional averaging method over 100 trials in subject I. This figure shows a bilateral beta rebound pattern with contralateral (left hemisphere) dominance.
with respect to the sensor array was determined by measur-ing magnetic signals from four head position indicator (HPI) coils placed on the scalp. Coil positions were identi-fied with a three-dimensional digitizer with respect to three predetermined landmarks (naison and bilateral preauricular points) on the scalp, and these data used to superimpose MEG source signals on individual MRI images obtained with a 3.0-T Bruker MedSpec S300 system (Bruker, Kalsrube, Germany). The anatomical image was acquired using a high-resolution T1-weighted, 3D gradient-echo pulse sequence (modified driven equilibrium Fourier
trans-form (MDEFT); TR/TE/TI ⫽ 88.1 ms/4.12 ms/650 ms,
128ⴱ128ⴱ128 matrix, FOV ⫽ 250 mm).
Empty room measurements were recorded for 3 min. Approximately 100 EOG-free trials of right index finger lifting were acquired and analyzed off-line. Since the focus was on beta-activities, the signals were further bandpass-filtered between 6 and 50 Hz (zero-phase, 10th-order, IIR Butterworth filter) to remove dc drifts and 60-Hz noise. The initial finger movement (movement onset; zero time) was registered with an optical switch (Taniguchi et al., 2000).
Electromyographic (EMG) activity from the extensor digi-torum communis (digitized at 1 KHz) was continuously recorded to monitor performance (see above). Each epoch comprised data points from⫺4 to 3 s relative to the move-ment onset (Salmelin et al., 1995; Salmelin and Hari, 1994a) and epochs were subjected to further single-trial ICA anal-ysis.
For SEF measurement, the right median nerve was elec-trically stimulated every 2 s with constant current pulses (0.3 ms in duration) exceeding the motor threshold. Approx-imately 100 EOG-free trials were acquired and digitized at 1 kHz and stored on magnetooptic disks for off-line analy-sis.
Data analysis
Independent component analysis of the single-trial MEG epoch
We take the advantages of sensitivity and localizing power of superficial sources by planar gradiometers (Rosell et al., 2001; Kajola et al., 1991). Each single-trial MEG
Fig. 6 (continued)
epoch contains m channels (m ⫽ 204, 102 pairs of
gradi-ometers) and n time points (usually m ⬍ n). The paired
gradiometer signals (⭸Bz/⭸x and ⭸Bz/⭸y) are arranged into
two (m/2)⫻ n submatrices B1and B2and concatenated into
an m ⫻ n matrix B. The ith rows (i ⱕ 102) of B1and B2
contain the measured gradiometer signals from the ith sen-sor location, and the jth column in B contains the measured data at the jth time point across all gradiometer channels. Mathematically, we can consider each row of B as samples generated from one random variable bi, i⫽ 1, 2,. . ., m. In
other words, matrix B is a realization of a random vector b
⫽ [b1b2. . .bm] T
.
The ICA techniques (Jung et al., 2001; Hyvarinen et al.,
2001) seek to find a p ⫻ m (p ⱕ m) matrix, W, which
converts the random vector b into another vector variable, s, consisting of p mutually independent random variables; thus; s p⫻1 ⫽
冤
s1 s2 · · · sp冥
⫽ W p⫻m b m⫻1 . (1)The mutual independence of si, for i⫽ 1, . . ., p, implies
that if P(si) represents the probability distribution of the ith
component, the joint probability distribution for all compo-nents can be factorized as
P共s1, s2, . . ., sP兲 ⫽ P共s1兲 P共s2兲 . . . P共sp兲. (2)
The ICA techniques use this assumption of mutual inde-pendence to find the unmixing matrix W.
All calculations in the present study were carried out using the FastICA algorithm which features high-speed cal-culation (cubic convergence) and does not require selection of step-size parameters or learning rate, unlike the gradient-based algorithm (Hyvarinen and Oja, 1997; Hyvarinen et al., 2001). The FastICA technique first removes means of row vectors in the B sample matrix such that each random variable bihas a zero mean, and then employs a whitening
process using principal component analysis. After whiten-ing, the covariance matrix of the whitened data becomes an identity matrix, and only the first p (pⱕ m) most significant principal components are preserved in the FastICA calcu-lation.
The next step is to look for a matrix that transforms the whitened data into a set of components as mutually inde-pendent as possible. Mutual information, as a measure of the independence of random variables, is used as the crite-rion for finding such a transformation. Mutual information can be expressed in terms of negentropy, an important measure of non-Gaussianity (Hyvarinen and Oja, 1997; Hyvarinen et al., 2001). Therefore, the problem of finding the independent components (s) and the transform matrix (W) can be translated into a search for linear combinations
of the whitened data that maximize the negentropy of the distributions of si, for i ⫽ 1, . . ., p.
After applying FastICA to the preprocessed single-trial MEG epochs, matrix B can be factored into a (mixing) matrix U and an (independent source) matrix S as
Bm⫻n⫽
冋
B1 B2册
⫽ Um⫻pSp⫻n ⫽冤
冋
u1,1 · · · u1,p · · · · · · um/ 2,1 · · · um/ 2,p册
冋
um/ 2⫹1,1 · · · um/ 2⫹1,p · · · · · · um,1 · · · um,p册
冥
m⫻p冤
sជ1 s ជ2 · · · s ជp冥
p⫻n , (3)in which each row sជiof matrix S僆ℜp⫻nrepresents samples
of an independent component (IC) si, for i⫽ 1, . . ., p and
U僆ℜm⫻pis the pseudoinverse of matrix W whose column
vectors represent the weight distribution values of the cor-responding ICs in S across all MEG gradiometer channels. In fact, matrix U is the “mixing matrix” that combines the
p ICs to reconstruct signal B. These temporal ICs can be
categorized into task-related ICs and task-unrelated ICs. Since the elicited brain activities or artifacts can be distrib-uted over multiple ICs, no one-to-one correspondence be-tween IC and source information is projected (Makeig et al., 1997). To facilitate the selection of task-related ICs, a tem-poral and spatial template pair was constructed prior to selection (see below). Spatial map xជj of the jth IC was
defined as the topographic display of all vector norms for weights of 102 gradiometer pairs in the jth column vector of
U, x ជj⫽ 关
冑
u1, j 2 ⫹ u 共m/ 2⫹1兲, j 2冑
u2, j 2 ⫹ u 共m/ 2⫹2兲, j 2 · · ·冑
um/ 2, j 2 ⫹ u m, j 2 兴T , (4)in which ui,jis the entry in the ith row and jth column of U
in Eq. (3). The spatial map is intended for component selection (see below).
Creation of a temporal template (VAMWtemplate) using amplitude modulation (envelope) of the MEG data
The recorded MEG signals at each gradiometer are fil-tered in the task-specific frequency band (Pfurtscheller and Lopes da Silva, 1999) and rectified by computing the AM waveform (envelope) using the amplitude modulation (AM) method (Clochon et al., 1996) as
m共t兲 ⫽
冑
MBP共t兲2⫹ H共MBP共t兲兲2, (5)in which MBP(t) is the bandpassed MEG signal and
H(MBP(t)) is its Hilbert transform. The task-specific fre-quency band is determined by the contrast between two 1-s amplitude spectra calculated over about 100 event-related EEG trials (Pfurtscheller and Lopes da Silva, 1999). One (serving as rest reference) is computed over the duration from 4 to 3 s preceding the onset of movement, and the
other (serving as reactive target) from 0.8 to 1.8 s after theonset of movement (see Figs. 1a and b). All beta-fre-quency components with significant modulation in terms of postmovement amplitude increase (above 95% confidence level; i.e., Z⬎ 3.09, P ⬍ 0.01) in the differential amplitude spectrum (see Fig. 1c) are taken as the task-specific fre-quency band for subsequent processing (Pfurtscheller and Berghold, 1989).
The vector norm of AM waveforms (VAMW) at each sensor site is computed using the square root of the AM
waveforms of each gradiometer pair; i.e., V (i, t) ⫽
公mx(i, t)2⫹ my(i, t)2, in which V(i, t) is the VAMW at the ith sensor location, and mx(i, t) and my(i, t) are the AM
waveforms in⭸Bz/⭸x and ⭸Bz/⭸y directions of the ith sensor
location. Event-related beta modulation is then computed as the difference in amplitude between the maximum ampli-tude of VAMW for each sensor site in the postmovement (0.8 to 1.8 s) interval and mean activity between⫺2.5 and
⫺2 s (see Fig. 2a) (Leocani et al., 1997). Beta rebound (BR)
is defined as the maximum amplitude of the computed event-related beta modulation from the subset of nine sensor sites in the vicinity of SMI (identified by SEF). The VAMWs of the BR calculation were averaged across the subjects (500 trials, 100 trials for each subject, five subjects pooled) to create the common temporal template, desig-nated VAMWtemplate(Fig. 2a).
Creation of a spatial template using topographical distribution of event-related beta modulation values
Individual spatial templates were first generated from the topographical distributions of event-related beta modulation values (see above). The five templates from the model generation group were then averaged to generate a common spatial template. In order to optimize conditions for spatial averaging, subjects’ heads were carefully positioned before actual measurements to keep head positioning and orienta-tion as similar as possible. Distances between head centers of the five subjects and the reference point (the origin of the MEG sensor array) in the horizontal plane were less than 4 mm, and angles between the vertical axis of the helmet and that of the head (the normal vector of the plane constituted by the three landmark points, i.e., nasion, and both preau-ricular points) remained within 5.5° (maximum deviation 1.5°) between subjects. Only the left half of the spatial map (unshaded in Fig. 2b) was used as the spatial template because this study focused on beta event-related activities in the hemisphere contralateral to the side of finger lifting; however, the other half can be generated analogously to extract activities in the ipsilateral hemisphere. Correlations among individual spatial templates ranged from 0.92 to 0.68. Respective correlations between the common spatial template and the individual spatial templates were 0.973, 0.811, 0.881, 0.904, and 0.915. These high correlation val-ues support the use of the spatial template in component selection for each individual’s magnetic signals.
Selection of pertinent independent components for the reconstruction of reactive beta activities
A spatial map (Eq. (4)) and corresponding VAMWs of each IC were generated for the selection of task-related ICs. Since the original signals may be decomposed into multiple ICs, the spectrum of each IC may vary from the one in the original signal due to the decomposition process. When settings for bandpass filtering for VAMW computation can-not be optimally determined using two-spectrum compari-son for the generation of a VAMWtemplate(Pfurtscheller and Lopes da Silva, 1999), three standard beta bands, 12–16, 16 –20, and 20 –24 Hz (Pfurtscheller, 1981), enclosing the event-related beta activities in motor task, were used to bandpass-filter (zero-phase, 10th-order, IIR Butterworth fil-ter) for each single-trial IC such that the three
frequency-laden resultant VAMWICs (the VAMWs bandpass-filtered
in three frequency bands of each IC) retained all task-related
information. These VAMWICs were subsequently used in
the selection of task-related ICs, which must fulfill the following dual criteria: (1) at least one of three
correspond-ing VAMWICs has a correlation with the VAMWtemplate
higher than 95% (Z⬎ 1.63, P ⬍ 0.05) among VAMWICs of
all the ICs for that single epoch, and (2) correlation between
the spatial map and spatial template is above 95% (Z ⬎
1.63, P ⬍ 0.05) for the spatial maps of all ICs. Data
processed via three-standard band filtering are not used in subsequent data reconstruction, but rather are used in con-junction with the dual-criteria only in the procedure “select-ing” the pertinent ICs. Unselected columns, i.e., task-unre-lated components, of mixing matrix U (Eq. (3)) are zeroed to produce a matrix Uˆ such that task-related rhythmic sig-nals are reconstructed by multiplying Uˆ and S (Fig. 3). The reconstructed data in each trial are then filtered within a trial-specific frequency band to extract reactive beta activ-ities.
Detection of task-laden trial-specific frequency band and extraction of reactive beta activities
The trial-specific frequency band detected in each trial is used to confine the reconstructed data within the most re-active beta band for further BR computation and source estimation. This frequency band is defined by the reactive beta band of the sensor site showing highest event-related beta modulation value (see Creation of temporal template) over the nine SMI vicinal sensor sites (identified by SEF) and is identified using the aforementioned two-spectrum procedure which has been suggested to be the best approach for the determination of reactive frequencies (Pfurtscheller and Lopes da Silva, 1999). Following data filtering with a trial-specific frequency band (zero-phase, 10th-order, IIR Butterworth filter), reactive beta activities in each single epoch can be extracted. The extracted reactive beta activi-ties are then subjected to source estimation and BR compu-tations.
Fig. 7. Smearing of MEG profile and decrease of BR magnitude due to latency jittering. (a) Raster plot of normalized VAMWrecon_max
Sas sorted by the latency measured between the time of peak beta rebound and the movement onset. Black dashed line indicates movement onset time. (b) Latency jittering resulting in a smearing of the MEG profile and a decrease of BR magnitude when more VAMWrecon_max
Fig. 9. Comparisons of magnetic fields and source locations preprocessed with ICA bandpass trial-specific (upper panel) and task-specific bandpass filtering (lower panel). (a) Neuromagnetic field maps. Data preprocessed with ICA trial-specific bandpass filter (15.57⫾ 3.21⬃22.17 ⫾ 3.3 Hz) gives a much less noisy neuromagnetic field pattern than that processed with the task-specific bandpass filtering method (15–21 Hz) (Pfurtscheller and Lopes da Silva, 1999). Black vertical lines in the tracings of the left column denote time points of the corresponding field maps in the right column. (b) Source localizations by ECD model. Only dipoles in postmovement interval-of-interest (IOI) with goodness-of-fit higher than 80% are accepted. The one with highest goodness-of-fit value out of each trial is rendered onto the subjects’ 3D MRI surfaces. The estimated source positions preprocessed by ICA bandpass filtering (upper panel) are (x, y, z)⫽ (⫺45 ⫾ 4.45, ⫺3.9 ⫾ 6.33, 80.7⫾ 3.63 mm; goodness-of-fit ⫽ 97.5 ⫾ 3.7%) in subject I (65 trials) and (x, y, z) ⫽ (⫺35.3 ⫾ 3.5, 5.7 ⫾ 6.02, 88.7 ⫾ 5.61 mm; goodness-of-fit ⫽ 96.9 ⫾ 3.7%) in subject VI (71 trials), whereas task-specific bandpass filtering (lower panel) yields (x, y, z)⫽ (⫺46.3 ⫾ 11.6, ⫺9.99 ⫾ 10.3, 84.51 ⫾ 6.7; goodness-of-fit
⫽ 89.7 ⫾ 3.4%) in subject I (65 trials) and (x, y, z) ⫽ (⫺31.8 ⫾ 8.13, 0.5 ⫾ 14.01, 87.9 ⫾ 12.54 mm; goodness-of-fit ⫽ 87.2 ⫾ 4.5%) in subject VI (71 trials),
respectively. The ICA trial-specific bandpass procedure yields better results in terms of more focused source locations and higher goodness-of-fit. x, y, and z denote the dipole location in the head coordinate system as anchored by the head position indicator (HPI) coils. The x axis passes through the preauricular points, pointing to the right; the positive y axis traverses the nasion and is normal to the x axis; the positive z axis points upward and normal to the xy plane.
Calculation of VAMWreconof reactive beta activities and single-trial epoch selection using a nonparametric sign test
Movement-related BR can be quantified from
single-epoch reactive beta activities and VAMWrecon(VAMW of
reconstructed data) for reactive beta activity at each sensor
site computed. The VAMWrecon of highest event-related
beta modulation (see Creation of temporal template) among the nine sensor sites vicinal to SMI is designated as VAMWrecon_max and is used in turn for single-trial epoch selection and BR computation, as the sensor site expressing VAMWrecon_maxdid not change throughout the experiment in our observations. A deterministic procedure, modified from Brovelli’s et al. (2002) approach, is used to select the significant trial. A nonparametric sign test is applied to the VAMWrecon_maxdesignated for BR calculation in each sin-gle trial by computing the Z score at each time point as
Z(t)⫽ (N⫹(t)⫺ (1/2)N)/((1/2)公N), in which N⫹(t) denotes the number of trials whose magnitudes are larger than the median value of their baseline activities at time point t, and
N the total number of trials. Time points with Z values
greater than 3.09 (P⬍ 0.01) are defined as the time interval-of-interest (IOI). After the determination of IOI for each subject, another sign test is then applied to find epochs showing significant increases in amplitude (Z⬎ 1.63, P ⬍ 0.05) using ZIOI(i)⫽ (NIOI⫹ (i)⫺ (1/2)NIOI)/((1/2) 公NIOI), in which ZIOI(i) is the Z value of the ith trial, NIOI⫹ (i) is the number of data points in postmovement IOI with values larger than the median of baseline activities of the ith trial, and NIOIis the total number of time points in postmovement IOI (Brovelli et al., 2002). An example of single-trial epoch selection is given in Fig. 4 (Subject I). The first trial in Fig. 4 with a ZIOIscore equal to⫺4.53 is marked as an insig-nificant epoch and eliminated from further analysis.
Source estimation of the reactive beta activities
Source estimation of the MEG reactive beta activities was done using equivalent current dipole (ECD) analysis and minimum current estimation (MCE) (Uutela et al., 1999; toolbox provided by Neuromag Ltd., Helsinki, Fin-land). A single dipole model was applied to explain the field every 1 ms, and only dipoles showing goodness-of-fit (Jensen and Vanni, 2002) values higher than 80% were used for data explanation. In MCE, the lattice constant of the triangular grid was 10 mm and locations closer than 30 mm to the center of the conductor were excluded from current estimates. Both analyses used a realistic head model for each subject. Template generation and single-trial data pro-cessing procedure are schematized in Figs. 5a and 5b, re-spectively. Epochs achieving significance in the increase of beta activities were chosen for subsequent BR calculation and dipole/source analysis.
Validation of coupled common spatial and temporal templates for single-trial analysis
Since there are inevitably differences in head size and variations in head positions inside the MEG scanner among
subjects, BR amplitude differences were compared using both individual spatial templates and the common spatial template. The use of a pair of common spatial and temporal templates for the extraction of individuals’ neuromagnetic single-trial signals was further validated on one additional subject.
Results
Due to space constraints, some results are illustrated using data from one subject. However, similar findings were observed in all subjects. Based on the known spatial loca-tion and temporal expression in terms of spatial and tem-poral templates, reactive beta activities were successfully extracted. Fig. 3a shows that IC 3 meets the dual criteria: (i) the correlation values between spatial map and spatial
tem-plate is 0.84 (rank ⫽ 97%, Z ⫽ 1.89, P ⫽ 0.03); (ii)
correlation values of 16 –20 Hz VAMWIC and 20 –24 Hz
VAMWICvs VAMWtemplateare 0.8 (rank⫽ 99%, Z ⫽ 3.08,
P⫽ 0.01) and 0.78 (rank ⫽ 97.8%, Z ⫽ 2.85, P ⫽ 0.022),
respectively. Fig. 3a illustrates that noise could also be identified and removed. IC 2 in Fig. 3a correlates highly (⫽0.88) in spatial distribution with the IC1 extracted from empty room measurements (Fig. 3b), and is therefore re-jected.
Fig. 6a depicts the single-trial VAMWrecons of subject I filtered within the trial-specific frequency band (Fig. 6b). The conventional AM method on the average of 100 epochs reveals a bilateral postmovement rebound pattern with con-tralateral (left hemisphere) dominance, whereas the current ICA-based single-trial analysis (one hemisphere template) yields only activation (one trial) in the left hemisphere (Figs. 6a and c).
Epoch acceptance rates were 84% (65/78), 89% (83/91), 71% (60/85), 73% (68/93), and 87% (76/87), respectively, for the model generation group and 81% (71/88) for the validation subject; the average for all six was 80.8%. The IOIs of significance were 0.76 –2.1 s, 0.66 –1.5 s, 0.8 –1.75 s, 0.46 –1.49 s, and 0.71–1.28 s for the five subjects in the model generation group, and 0.88 –1.67 s for the validation subject. Averaged magnitude of BR was calculated from the reconstructed data on trials that survived the epoch-selec-tion procedure. The BR amplitudes computed from individ-ual spatial templates were 20.9⫾ 7.1 (mean ⫾ SD), 18.1 ⫾ 10.3, 16.2⫾ 6.2, 23.2 ⫾ 10.89, and 6.2 ⫾ 2.7 for the first five subjects, respectively, and 27.6 ⫾ 11.1 fT/cm for the sixth subject (Table 1). Using the common spatial template, BR amplitudes were 21.1⫾ 7.97, 19.02 ⫾ 9.7, 15.5 ⫾ 5.3, 19.75⫾ 8.75, 5.91 ⫾ 3.2, and 27.1 ⫾ 10.2 fT/cm, respec-tively (Table 1). There was no significant difference be-tween the results obtained with two approaches (P⫽ 0.88; unpaired two-tailed t test). BR amplitudes obtained with the conventional method of averaging on 100 trials were 18.2, 7.254, 12.92, 16.4, 2.9, and 23.12 fT/cm, respectively. Means for single-trial ICA-derived BRs, using either
indi-vidual or common spatial templates, were significantly higher than those obtained using the conventional method of
averaging (P⬍ 0.005; matched-pair Wilcoxon test; Table
1). The comparisons of BR amplitude and task-specific frequency band between ICA-based single-trial and conven-tional methods are given in Table 1.
The ICA-based single-trial approach shows remarkable latency jittering and intertrial variability throughout the whole measurement process. Both factors can result in at-tenuation and smearing of averaged movement-related MEG responses. Fig. 7a shows the raster plot of 65 normal-ized single-trial VAMWrecon_maxs which survived the selec-tion procedure for subject I, sorted by VAMWrecon_maxpeak latency as indexed to movement onset. The mean latency of peak beta rebound for the 65 trials was 1.41⫾ 0.43 s (mean
⫾ SD). With more epochs (random selection) averaged as
with the conventional method of averaging, the averaged BR was attenuated (25.3, 24.6, 22.3, 21.5, 20.3, and 21.1 fT/cm for 1, 15, 25, 35, 45, and 65 trials averaged,
respec-tively; values taken from the averaged VAMWrecon_maxs
using common spatial template) and the time–activity plots were smeared (Fig. 7b).
Source estimation using ECD and MCE both showed a cluster of current sources centered (mean coordinates) in the anterior bank of the central sulcus (see Figs. 8e and f) on data points around the rebound peak of extracted reactive beta activities (see Fig. 8d, time interval between 1202 and 1302 ms of one single epoch of subject I). The ECD-located dipoles oscillate and span a sector. Furthermore, the center of MCE-estimated current sources (yellow dots) lies less than 2 mm from the center of ECD-estimated dipoles (red dots) (see Fig. 8f). These results cross-verify the validity of the ICA-based single-trial method.
Discussion
The movement-related oscillatory modulations (ERD/ ERS of alpha, beta, and gamma) have been reported to be spatially extended (Babiloni et al., 1999; Crone et al., 1998a; Crone et al., 1998b; Leocani et al., 1997; Neuper and
Pfurtscheller, 2001; Salmelin and Hari, 1994a; Taniguchi et al., 2000; van Burik et al., 1998). Source localizations using conventional filtering have also been reported to disperse among several regions (Salmelin and Hari, 1994a). How-ever, our results strongly indicate that proper treatment when trial-by-trial dynamics can be accounted for yields clustered localizations congruent to neuroanatomical repre-sentations.
The present ICA-based spatiotemporal approach for sin-gle-trial analysis study is dedicated to the extraction of neuromagnetic measurements of event-related beta oscilla-tory activities. One distinct feature of the current ICA-based method compared with other single-trial approaches (Guger et al., 2000; Ioannides et al., 1993; Jung et al., 2001) is the simultaneous use of a spatial template and a temporal tem-plate for component selection. The spatial temtem-plate provides a priori spatial information for brain signals, while the temporal template contains temporal characteristics of event-related responses. Using the paired criteria for com-ponent selection, identification specificity of task-related components for signal reconstruction is significantly im-proved. As shown in Fig. 3c, the inclusion of IC 9 with high spatial correlation (correlation value⫽ 0.61, rank ⫽ 95.2%,
Z ⫽ 1.67, P ⫽ 0.048) but devoid of temporal congruence
(correlation value⫽ 0.28, rank ⫽ 13%, Z ⫽ 0.34, P ⫽ 0.87) causes beta BR to deteriorate from 28.9 fT/cm (red curve) to 18.6 fT/cm (blue curve). The ICA-preprocessed data set yields cleaner field maps (Fig. 9a), which result in circum-scribed localizations (Fig. 9b; cf., Salmelin and Hari, 1994a).
Significantly, the current method also makes possible the analysis of the reactive frequency band for every single trial once task-related rhythmic activities are extracted. The con-ventional method discounts this subtle but potentially im-portant information. Notwithstanding, the idea of using a fixed window for signal filtering is neurophysiologically not optimal. We emphasize the precise identification of reactive trial-specific frequencies for BR calculation, since task-related frequency modulation might exist in one or multiple bands (Pfurtscheller and Lopes da Silva, 1999). The three-standard frequency band procedure is used for generation of
Table 1
Comparison of BR amplitude and specific frequency bands for ICA-based single-trial and conventional methods
Subject index ICA based single-trial method Conventional AM method
BR amplitude (fT/cm) Trial-specific frequency band (Hz) BR amplitude (fT/cm) Task-specific frequency band (Hz) Individual spatial template Common spatial template Individual spatial template Common spatial template I 20.9⫾ 7.1 21.1⫾ 7.97 16.67⫾ 2.77–21.22 ⫾ 2.44 15.57⫾ 3.21–22.17 ⫾ 3.3 18.2 15–21 II 18.1⫾ 10.3 19.02⫾ 9.7 18.04⫾ 2.62–22.18 ⫾ 3.12 17.92⫾ 2.3–21.9 ⫾ 2.72 7.25 17–20 III 16.2⫾ 6.2 15.5⫾ 5.3 16.2⫾ 1.89–20.49 ⫾ 2.3 16.8⫾ 2.3–20.91 ⫾ 2.22 12.92 15–19 IV 23.2⫾ 10.89 19.75⫾ 8.75 16.1⫾ 2.37–20.7 ⫾ 3.08 15.5⫾ 3.3–19.2 ⫾ 2.77 16.4 14–17 V 6.2⫾ 2.7 5.91⫾ 3.2 17.31⫾ 3.23–20.77 ⫾ 3.67 16.8⫾ 3.1–21.2 ⫾ 2.9 2.9 17–20 VI (validation) 27.6⫾ 11.1 27.1⫾ 10.2 16.32⫾ 2.83–19.94 ⫾ 2.68 16.81⫾ 2.72–20.14 ⫾ 3.1 23.12 16–20
Fig. 8. Overlay of extracted reactive beta activities on MR image. (a) Spatial map reconstructed using xជrecon⫽ ⌺jk⫽1xជj(see Materials and methods). (b, c)
Different views of superposition of the isocontour spatial map on the segmented MRI brain. (d) Representative trace of reconstructed reactive beta activities in the vicinity of SMI. (e) Upper panels are isocontour maps of reconstructed neuromagnetic signals at 1202 mspost movement. Lower panels show that all dipoles (from 1202 to 1302 ms after movement onset as box-framed in (d)) are located in the primary motor area and oscillate accordingly. (f) The center of the MCE-estimated current sources (yellow dot) overlays the source location determined using the equivalent current dipole method (ECD) (red dot). Upper left panel: coronal view; upper right panel: sagittal view; lower left panel: axial view; lower right panel: distribution of MCE estimated current sources.
VAMWICs to recover all possible task-related information and is followed by a two short-time spectra comparison procedure (Pfurtscheller and Lopes da Silva, 1999) for the identification of the optimal reactive trial-specific frequency band in the reconstructed epochs. The present approach not only extracts the specific reactive frequencies but also re-tains phase information on a by-trial basis. The trial-specific frequency band of postmovement beta modulation anchors mainly (⬃85% of all trials) in the lower beta band (⬃16–20 Hz) and less frequently (⬃15%) in the higher beta band (20 –24 Hz). Great variation of BR values is also seen, as reflected in large SD (Table 1). The revealed trial-by-trial dynamics provide a possibility for future profound study of subtle brain dynamics.
It is noteworthy that not all the data reconstructed from the selected ICs survive the statistical threshold. We have carefully monitored online and thoroughly checked offline
the EMG measurements in terms of EMG onset (P⫽ 0.61,
unpaired two-tailed t test), termination (P⫽ 0.53, unpaired
two-tailed t test) and the EMG duration (P ⫽ 0.573,
un-paired two-tailed t test) during finger lifting between sig-nificant and insigsig-nificant trials as indexed to the movement
registration by the optic pad (Abbink et al., 1998). The data indicate an absence of prominent behavioral difference commensurate to the differential neuromagnetic responses. Some epochs with a fluctuating baseline, e.g., non-task-related spontaneous bursts of beta oscillatory activities, may manifest high baseline activity, which in turn results in a decrease in BR readout leading to exclusion after statistical manipulation (Fig. 4). It has been suggested that baseline spontaneous activities may carry important information rel-evant to attention level, wakefulness, task difficulty, etc. (Buser and Rougeul-Buser, 1999; Sterman, 1999). The jit-tering of the neuromagnetic beta ERS is likewise interesting and may be also physiological. A zero-phase Butterworth filter was used to bandpass-filter the raw data. The symmet-ric property of the zero-phase filter means that processed signals have precisely zero phase distortion and therefore no time shift of peak beta rebound was introduced. Hence, fluctuations of significance level and the jittering of central processing despite similar behavioral performance may be ascribed to the subject’s variant cognitive states or the degree of training (Buser and Rougeul-Buser, 1999; Ster-man, 1999; Flotzinger et al., 1992; Wolpaw and McFarland,
Fig. 8 (continued)
1994; Bastiaansen et al., 1999, 2001; Earle, 1988; Haig et al., 1995; Hoffman et al., 1991; Yabe et al., 1993). The exploration of underlying mechanisms mandates more me-ticulous designs in the future. Using the conventional method of averaging, certain diseases, such as Parkinson’s and Unverricht–Lundborg myoclonic epilepsy, have been observed to show attenuated, prolonged, or abolished ERS responses (Silen et al., 2000; Tamas et al., 2003). Such cases can be further examined using the current ICA-based single-trial method for the time course and single-trial-by-single-trial dynamics to disclose hitherto unexplored mechanisms underlying these phenomena.
A concern with any data-driven method is that prominent artifacts or noise can be intermingled with task-specific information (Ermer et al., 2000; Lins et al., 1993a, 1993b). However, previous ICA reports (Makeig et al., 1997; McKeown and Radtke, 2001) indicate that brain rhythmic signals generated from different sources usually have their own oscillatory frequencies with distinct phases and are located in specific brain regions with patterns that are dis-tinct from artifacts or noise (see also Fig. 3). This endorses the feasibility of using ICA to separate targeted rhythmic signals from irrelevant ones. The high epoch-acceptance
rate (⬃80%) can be attributed to an improved SNR
com-pared to other studies on single-trial approaches to sensori-motor oscillatory activities (Brovelli et al., 2002; Wolpaw and McFarland, 1994). For instance, the spatial map of IC2 in Fig. 2a correlates highly (0.88) with the spatial map of IC1 from empty room measurement as shown in Fig. 3b; this suggests that the neuromagnetic signal IC2, deselected for subsequent processing, can be accounted for by back-ground noise in the shielding room. IC11 in Fig. 3a has a stationary cycle around 1.2 Hz, and its spatial map has higher weights at the outer rim of the MEG sensor array, which suggests a plausible connection with cardiac cycles. It was also observed (Fig. 3) that rhythmic activities in left and right SMIs as well as the occipital areas could be extracted into separate ICs that can be reminiscent of vari-ous mechanisms and time courses of different brain oscil-latory activities (Pfurtscheller and Lopes da Silva, 1999; Pfurtscheller et al., 1997, 1998b; Stancak and Pfurtscheller, 1996a, 1996b; Andrew and Pfurtscheller, 1999).
Since most task irrelevant signals, e.g., internal and ex-ternal noises, can be removed by proper deselection of ICs, it is possible to reconstitute the representative spatial map of all contributing ICs using xជrecon⫽ ⌺j⫽1
k x
ជj, in which xជreconis the reconstructed spatial map, k is the number of selected ICs, and xជiis the spatial map of the ith selected IC in Fig. 4
(Fig. 8a). This spatial map of reconstructed signals, which is a topographical distribution of weighting factors on the sensor array, can be overlaid with the segmented MRI brain (Figs. 8b and c; ASA program, ANT Software, Dutch). The highest weight is shown to project over the SMI area, which demonstrates that the high SNR of the ICA-extracted rhyth-mic activities of each trial has made possible the use on single-trial data of source estimation methods that require
high SNR on input data for processing, e.g., the ECD technique, MCE, and MNE (Delorme et al., 2001, 2002; Jung et al., 2001; Makeig et al., 1997; McKeown and Radke, 2001). Conventionally, these estimation methods exploit averaged data out of a large amount of trials.
Another reason why the intricate phase-unlocked signal can be preserved is the fact that no averaging procedure is needed; such a procedure would otherwise inherently distort the embedded information. Accordingly, as shown in Fig. 8d, source modeling with a moving dipole on a millisecond-by-millisecond basis on the reconstructed oscillatory beta signals during the rebound period (Brovelli et al., 2002) of a single-trial epoch results in a focused clustering of dipole foci at the precentral area, i.e., the primary motor cortex (Fig. 8e). Fig. 8 shows the result of MCE modeling (Uutela et al., 1999), where the center of MCE-estimated current sources (yellow dot) is very close (⬍2 mm distance) to the dipole location as estimated using the ECD approach (red dot).
It can be argued that one can first localize the generator area and then build a spatial filter for extracting single-trial data so that the subsequent analysis can be conducted on the source level instead of the sensor level. One premise and a justification of using a source-area-generated spatial filter are that the source area can be precisely localized for the generation of a spatial filter (Tesche et al., 1995). The very first step is to filter the signals to obtain a presupposed reactive frequency band. However, using conventional sim-ple filtering techniques, ambient noise with ⬃20-Hz com-ponents cannot be optimally removed, and this will cause localization uncertainty for the probed sources (Fig. 9). However, ICA preprocessing decomposes the compound neuromagnetic signals into various independent task-related and task-unrelated/noise components so that⬃20-Hz activ-ities not related to the a priori spatiotemporal profile will not confound the selected ones. Furthermore, our ICA-based method differs from other spatial filtering techniques, e.g., signal space projection (SSP) which is a fixed spatial filter for signal extraction (Tesche et al., 1995). The ICA-based method blindly decomposes the MEG epochs (B) into a spatially distributed map (U) multiplied by temporal signals (S), i.e., B ⫽ U · S, on the basis of independency among sources (Vigario and Oja, 2000), whereas SSP mandates a predefined spatial filter (U_sf) for recovering signals (S), i.e.,
S ⫽ U_sf⫹ · B, where ⫹ denotes pseudo inverse, based on orthogonal projection. When ambient noise and the spatial filter are not mutually orthogonal, the SSP has difficulty in resolving the two. Subsequent application of ICA following SSP does not ensure finer signal extraction or further noise removal since the data recovered from SSP are already linear mixtures of components of a predefined spatial filter, which is a constraint drag on the optimal performance of ICA designed for blind decomposition.
Left and right sensorimotor rhythms can be decomposed into two distinct ICs (IC 3 and IC 5 in Fig. 3), respectively, implying possible independent modulatory mechanisms
be-tween the two hemispheres. This view is corroborated by an event-related coherence study (Andrew and Pfurtscheller, 1999) that reports a lack of interhemispheric coherence in human postmovement beta activities. Movement-related beta oscillatory activities of the right hemisphere can be extracted in the same way using spatial and temporal tem-plates for right sensorimotor rhythm. The source locations for extracted right hemispheric beta activities were mainly in the right premotor area (data not shown), which agrees with previous studies (Brovelli et al., 2002; Ilmoniemi, 1991). Event-related beta activities in SMA and posterior parietal cortical areas (Brovelli et al., 2002; Joliot et al., 1999) are not observed in our data, possibly due to the fact that the contributing sources here are radial in orientation and thus could not be optimally detected by MEG (Salmelin and Hari, 1994b).
The agreement between the values of BR amplitude obtained with the common spatial/temporal templates and the individually generated ones (Table 1) promises a flexi-bility in both experimental design and analytical strategy. The proposed ICA-based spatiotemporal approach for sin-gle-trial analysis can also be applied on fewer trials (Fig. 7b), which is a great advantage over conventional methods. Given meticulous head positioning (see Materials and meth-ods), common spatial and temporal templates can be used to extract pertinent movement-related neuromagnetic signals from subjects, which may shorten the overall time needed to run an experiment. We have no preference for the use of a grand-averaged template over individual ones. On the con-trary, the use of an individual template is suggested for any profound individual-based ERD/ERS study. However, the feasibility of using a grand averaged template provides an effective alternative in cases where lengthy procedures can-not be endured by the participants. This is particularly true for clinical settings where patients have attention problems or are incapable of sustaining-long experiments so that individual templates cannot be optimally obtained. Never-theless, caution should be exercised when applying the current ICA-based single-trial method for clinical studies. For patients whose heads cannot be properly positioned in the center of the MEG helmet, the use of a common spatial template may fail, making a customized individual spatial template mandatory for IC selection. For patients whose motor performance deviates significantly from normal, e.g., victims of motor stroke or severe movement disorders, the use of the common temporal template might not be justified since the time courses of event-related brain activities may be significantly altered due to primary deficit or secondary plasticity. Accordingly, in such situations, an individual spatial template can be applied without a temporal template as an aid to component selection. Our future investigations will combine the current dual-template approach with a source estimation method so that a spatial filter of better precision and higher dimensions can be designed, which will make possible sophisticated analysis on the source level
instead of the sensor level, eliminating the positioning prob-lem.
Conclusions
The present novel ICA-based spatiotemporal approach for single-trial analysis features a paired-template matching for stringent component selection. The spatial template pro-vides a priori spatial information for targeted brain signals while the temporal template contains temporal characteris-tics of event-related responses. Our method takes into ac-count the subtle trial-by-trial dynamics. The method prom-ises not only a high extraction rate of postmovement beta synchronization but also better localization of the corre-sponding sources. Various source modeling methods com-manding high SNR can now be applied to single-trial data extracted using the ICA–spatiotemporal procedure. The em-bodied common template approach permits an effective alternative in cases where lengthy procedures cannot be endured by the participants or in clinical settings where patients have attention problems or are incapable of sustain-ing long experiments.
Acknowledgments
This study was funded by the Taipei Veterans General Hospital (90400, 90443, 91361, 91380, 923721, 92348, VTY-929112), National Science Council (902314B075124, 902314B075115, 902511S010001, 912314B075069, 912314B075069), and the Ministry of Education (89BFA221401) of Taiwan. We thank Mr. Chi-Cher Chou for MRI assistance and Dr. David M. Niddam for previous input to the work.
References
Abbink, J.H., van der Bilt, A., van der Glas, H.W., 1998. Detection of onset and termination of muscle activity in surface electromyograms. J. Oral Rehabil. 25, 265–269.
Andrew, C., Pfurtscheller, G., 1995. Event-related coherence during finger movement: a pilot study. Biomedizinische Technik 40, 326 –332. Andrew, C., Pfurtscheller, G., 1999. Lack of bilateral coherence of
post-movement central beta oscillations in the human electroencephalogra-phy. Neurosci. Lett. 273, 82– 89.
Babiloni, C., Carducci, F., Cincotti, F., Rossini, P.M., Neuper, C., Pfurtscheller, G., Babiloni, F., 1999. Human movement-related poten-tials vs desynchronization of EEG alpha rhythm: a high-resolution EEG study. NeuroImage 10, 658 – 665.
Bastiaansen, M.C.M., Bocker, K.B.E., Brunia, C.H.M., 2001. Event-re-lated desynchronization during anticipatory attention for an upcoming stimulus: a comparative EEG/MEG study. Clin. Neurophysiol. 112, 393– 403.
Bastiaansen, M.C.M., Bocker, K.B.E., Cluitmans, P.J.M., Brunia, C.H.M., 1999. Event-related desynchronization related to the anticipation of a stimulus providing knowledge of results. Clin. Neurophysiol. 110, 250 –260.
Brovelli, A., Battaglini, P.P., Naranjo, J.R., Budai, R., 2002. Medium-range oscillatory network and the 20-Hz sensorimotor induced poten-tial. NeuroImage 16, 130 –141.
Buser, P., Rougeul-Buser, A., 1999. EEG synchronization in cat, monkey and human during attentive states: a brief survey, in: Pfurtscheller, G., Lopes da Silva, F.H. (Eds.), Event-related desynchronization. Hand-book of electroencephalography and clinical neurophysiology, Elsevier, Amsterdam, pp. 13–32.
Clochon, P., Fontbonne, J.M., Etevenon, P., 1996. A new method for quantifying EEG event-related desynchronization: amplitude envelope analysis. Electroencephalogr. Clin. Neurophysiol. 98, 126 –129. Crone, N.E., Miglioretti, D.L., Gordon, B., Lesser, R.P., 1998. Functional
mapping of human sensorimotor cortex with electrocorticographic spectral analysis. II. Event-related synchronization in the gamma band. Brain 121, 2301–2315.
Crone, N.E., Miglioretti, D.L., Gordon, B., Sieracki, J.M., Wilson, M.T., Uematsu, S., Lesser, R.P., 1998. Functional mapping of human senso-rimotor cortex with electrocorticographic spectral analysis. I. Alpha and beta event-related desynchronization. Brain 121, 2271–2299. Delorme, A., Makeig, S., Fabre-Thorpe, M., Sejnowske, T., 2001. From
single-trial EEG to Brain Area dynamics. Neurocomputing 44, 1057– 1064.
Delorme, A., Makeig, S., 2003. EEG changes accompanying learned reg-ulation of 12-Hz EEG activity. IEEE Trans. Neural Syst. Rehab. Engin. 2, 133–136.
Duann, J.R., Jung, T.P., Kuo, W.J., Yeh, T.C., Makeig, S., Hsieh, J.C., Sejnowski, T.J., 2002. Single-trail variability in event-related BOLD signals. NeuroImage 15, 823– 835.
Earle, J.B., 1988. Task difficulty and EEG alpha asymmetry: an amplitude and frequency analysis. Neuropsychobiology 20, 95–112.
Ermer, J.J., Mosher, J.C., Huang, M., Leahy, R.M., 2000. Paired MEG data set source localization using recursively applied and projected (RAP) MUSIC. IEEE Trans. Biomed. Engin. 47, 1248 –1260.
Florian, G., Pfurtscheller, G., 1995. Dynamic spectral analysis of event-related EEG data. Electroencephalogr. Clin. Neurophysiol. 95, 393– 396.
Flotzinger, D., Kalcher, J., Pfurtscheller, G., 1992. EEG classification by learning vector quantization. Biomedizinische Technik 37, 303–309. Guger, C., Ramoser, H., Pfurtscheller, G., 2000. Real-time EEG analysis
with subject-specific spatial patterns for a brain-computer interface (BCI). IEEE Trans. Rehabil. Engin. 8, 447– 456.
Haig, A.R., Gordon, E., Rogers, G., Anderson, J., 1995. Classification of single-trial ERP sub-types: application of globally optimal vector quan-tization using simulated annealing. Electroencephalogr. Clin. Neuro-physiol. 94, 288 –297.
Hoffman, R.E., Buchsbaum, M.S., Escobar, M.D., Makuch, R.W., Nuechterlein, K.H., Guich, S.M., 1991. EEG coherence of prefrontal areas in normal and schizophrenic males during perceptual activation. J. Neuropsychiatry Clin. Neurosci. 3, 169 –175.
Hyvarinen, A., Karhunen, J., Oja, E., 2001. Independent Component Anal-ysis. Wiley, New York.
Hyvarinen, A., Oja, E., 1997. A fast fixed-point algorithm for independent component analysis. Neural Comput. 9 (7), 1483–1492.
Ilmoniemi, R.J., 1991. Estimates of neuronal current distributions. Acta Otolaryngol. 49 (Suppl), 80 – 87.
Ioannides, A.A., Singh, K.D., Hasson, R., Baumann, S.B., Rogers, R.L., Guinto, F.C., Papanicolaou, A.C., 1993. Comparison of current dipole and magnetic field tomography analyses of the cortical response to auditory stimuli. Brain Topogr. 6, 27–34.
Jensen, O., Vanni, S., 2002. A new method to identify multiple sources of oscillatory activity from magnetoencephalographic data. NeuroImage 15, 568 –574.
Joliot, M., Papathanassiou, D., Mellet, E., Quinton, O., Mazoyer, N., Courtheoux, P., Mazoyer, B., 1999. FMRI and PET of self-paced finger movement: comparison of intersubject stereotaxic averaged data. Neu-roImage 10, 430 – 447.
Jung, T.P., Makeig, S., Westerfield, M., Townsend, J., Courchesne, E., Sejnowski, J.T., 2001. Analysis and visualization of single-trial event-related potentials. Human Brain Mapp. 14, 166 –185.
Kajola, M., Ahonen, A., Hamalainen, M.S., Knuutila, J., Lounasmaa, O.V., Simola, J., Vilkman, V., 1991. Development of multichannel neuro-magnetic instrumentation in Finland. Clin. Physics Physiologic. Meas. 12, 39 – 44.
Kao, Y.H., Guo, W.Y., Wu, Y.T., Liu, K.C., Chai, W.Y., Lin, C.Y., Hwang, Y.H., Liou, A.J.K., Cheng, H.C., Yeh, T.C., Hsieh, J.C., Teng, M.M.H., 2003. Hemodynamic segmentation of MR brain perfusion images using independent component, thresholding and Bayesian esti-mation. Magn. Reson. Med. 49, 885– 894.
Kalcher, J., Pfurtscheller, G., 1995. Discrimination between phase-locked and non-phase-locked event-related EEG activity. Electroencephalogr. Clin. Neurophysiol. 94, 381–384.
Klimesch, W., Russegger, H., Doppelmayr, M., Pachinger, T., 1998. A method for the calculation of induced band power: implications for the significance of brain oscillation. Electroencephalogr. Clin. Neuro-physiol. 108, 123–130.
Leocani, L., Manganotti, C.T., Zhuang, P., Hallett, M., 1997. Event-related coherence and event-related desynchronization/synchronization in the 10 Hz and 20 Hz EEG during self-paced movements. Electroencepha-logr. Clin. Neurophysiol. 103, 199 –206.
Lins, O.G., Picton, T.W., Berg, P., Scherg, M., 1993a. Ocular artifacts in EEG and event-reated potentials: I. Scalp topography. Brain Topogr. 6, 51– 63.
Lins, O.G., Picton, T.W., Berg, P., Scherg, M., 1993b. Ocular artifacts in recording EEGs and event-related potentials: II. Source dipoles and source components. Brain Topogr. 6, 65–78.
Lopes da Silva, F.H., 1991. Neural mechanisms underlying brain waves: from neural membranes to networks. Electroencephalogr. Clin. Neuro-physiol. 79, 81–93.
Lopes da Silva, F.H., 1996. The generation of electric and magnetic signals of the brain by local networks. Comp. Human Physiol. 1, 509 –531. Makeig, S., Jung, T.P., Bell, A.J., Ghahremani, D., Sejnowski, T., 1997.
Blind separation of auditory event-related brain responses into inde-pendent components. Proc. Natl. Acad. Sci. USA 94, 10979 –10984. McKeown, M.J., Makeig, S., Brown, G.G., Jung, T.P., Kindermann, S.S.,
Bell, A.J., Sejnowski, T.J., 1998. Analysis of fMRI data by blind separation into independent spatial components. Human Brain Mapp. 6, 160 –188.
McKeown, M., Radtke, R., 2001. Phasic and tonic coupling between EEG and EMG demonstrated with independent component analysis. J. Clin. Neurophysiol. 18, 45–57.
Mosher, J.C., Lewis, P.S., Leahy, R.M., 1992. Multiple dipole modeling and localization from spatio-temporal MEG data. IEEE Trans. Biomed. Engin. 39, 541–557.
Muller-Gerking, J., Pfurtscheller, G., Flyvbjerg, H., 1999. Designing op-timal spatial filters for single-trial EEG classification in a movement task. Clin. Neurophysiol. 110, 787–798.
Neuper, C., Pfurtscheller, G., 2001. Event-related dynamics of cortical rhythms: frequency-specfic features and functional correlates. Int. J. Psychophysiol. 43, 41–58.
Pfurtscheller, G., 1981. Central beta rhythm during sensorimotor activities in man. Electroencephalogr. Clin. Neurophysiol. 51, 253–264. Pfurtscheller, G., Aranibar, A., 1977. Event-related cortical
desynchroni-zation detected by power measurements of scalp EEG. Electroencepha-logr. Clin. Neurophysiol. 42, 817– 826.
Pfurtscheller, G., Aranibar, A., 1979. Evaluation of event-related desyn-chronization (ERD) preceding and following voluntary self-paced movement. Electroencephalogr. Clin. Neurophysiol. 46, 138 –146. Pfurtscheller, G., Berghold, A., 1989. Patterns of cortical activation during
planning of voluntary movement. Electroencephalogr. Clin. Neuro-physiol. 72, 250 –258.
Pfurtscheller, G., Lopes da Silva, F.H., 1999. Event-related EEG/MEG synchronization and desynchronization: basic principles. [Review]. Clin. Neurophysiol. 110, 1842–1857.
Pfurtscheller, G., Neuper, C., Flotzinger, D., Pregenzer, M., 1997. EEG-based discrimination between imagination of right and left hand move-ment. Electroencephalogr. Clin. Neurophysiol. 103, 642– 651. Pfurtscheller, G., Pichler-Zalaudek, K., Ortmayr, B., Diez, J., Reisecker, F.,
1998a. Postmovement beta synchronization in patients with Parkin-son’s disease. J. Clin. Neurophysiol. 15, 243–250.
Pfurtscheller, G., Pregenzer, M., Neuper, C., 1994. Visualization of sen-sorimotor areas involved in preparation for hand movement based on classification of alpha and central beta rhythms in single EEG trials in man. Neurosci. Lett. 181, 43– 46.
Pfurtscheller, G., Stancak Jr., A., Neuper, C., 1996. Post-movement beta synchronization: a correlate of an idling motor area? Electroencepha-logr. Clin. Neurophysiol. 98, 281–293.
Pfurtscheller, G., Zalaudek, K., Neuper, C., 1998b. Event-related beta synchronization after wrist, finger and thumb movement. Electroen-cephalogr. Clin. Neurophysiol. 109, 154 –160.
Rosell, J., Casanas, R., Scharfetter, H., 2001. Sensitivity maps and sys-temrequirements for magnetic induction tomography using a plannar gradiometer. Physiologic. Meas. 22, 121–130.
Salmelin, R., Hamalainen, M., Kajola, M., Hari, R., 1995. Functional segregation of movement-related rhythmic activity in the human brain. NeuroImage 2, 237–243.
Salmelin, R., Hari, R., 1994a. Characterization of spontaneous MEG rhythms in healthy adults. Electroencephalogr. Clin. Neurophysiol. 91, 237–248.
Salmelin, R., Hari, R., 1994b. Spatiotemporal characteristics of sensori-motor neuromagnetic rhythms related to thumb movement. Neuro-science 60, 537–550.
Silen, T., Forss, N., Jensen, O., Hari, R., 2000. Abnormal reactivity of the
⬃20-Hz motor cortex rhythm in Unverricht Lundborg type progressive
myoclonus epilepsy. NeuroImage 12, 707–712.
Stancak Jr., A., Pfurtscheller, G., 1996a. The effects of handedness and type of movement on the contralateral preponderance of mu-rhythm desynchronization. Electroencephalogr. Clin. Neurophysiol. 99, 174 – 182.
Stancak Jr., A., Pfurtscheller, G., 1996b. Mu-rhythm changes in brisk and slow self-paced finger movements. Neuroreport 7, 1161–1164. Sterman, M.B., 1999. Event-related EEG response correlates of task
dif-ficulty, sleep deprivation and sensory distraction, in: Pfurtscheller, G., Lopes da Silva, F.H. (Eds.), Event-Related desynchronization. Hand-book of electroencephalography and clinical neurophysiology, Elsevier, Amsterdam, pp. 233–242.
Tamas, G., Szirmai, I., Palvolgyi, L., Takats, A., Kamondi, A., 2003. Impairment of post-movement beta synchronization in parkinson’s disease is related to laterality of tremor. Clin. Neurophysiol. 114, 614 – 623.
Tang, A.C., Pearlmutter, B.A., Malaszenko, N.A., Phung, D.B., 2002. Independent components of magnetoencephalography: single-trial re-sponse onset times. NeuroImage 17, 1773–1789.
Taniguchi, M., Kato, A., Fujita, N., Hirata, M., Tanaka, H., Kihara, T., Ninomiya, H., Hirabuki, N., Nakamura, H., Robinson, S.E., Cheyne, D., Yoshimine, T., 2000. Movement-related desynchronization of the cerebral cortex studied with spatially filtered magnetoencephalograpy. NeuroImage 12, 298 –306.
Tesche, C.D., Unsitalo, M.A., Ilmoniemi, R.J., Huotilainen, M., Kajola, M., Salonen, O., 1995. Signal-space projections of MEG data charac-terize both distributed and well-localized neuroal sources. Electroen-cephalogr. Clin. Neurophysiol. 95, 189 –200.
Uutela, K., Hamalainen, M., Somersalo, E., 1999. Visualization of mag-netoencephalography data using minimum current estimates. NeuroIm-age 10, 173–180.
van Burik, M., Pfurtscheller, G., 1999. Functional imaging of postmove-ment beta event-related synchronization. J. Clin. Neurophysiol. 16, 383–390.
Vigario, R., Oja, E., 2000. Independence: a new criterion for the analysis of the electromagnetic fields in the global brain? Neural Netw. 13, 891–907.
Wolpaw, J.R., McFarland, D.J., 1994. Multichannel EEG-based brain-computer communication. Electroencephalogr. Clin. Neurophysiol. 90, 444 – 449.
Wu, Y.T., Lee, P.L., Chen, L.F., Yeh, T.C., Hsieh, J.C., 2002. Single-trial quantification of imagery beta-band Mu rhythm in finger lifting task using independent component analysis (ICA), in: Proceedings, BioMag 13thInternational Conference on Biomagnetism, pp. 1045–1047. Wu, Y.T., Lee, P.L., Chen, L.F., Yeh, T.C., Hsieh, J.C., 2003.
Quantifi-cation of movement-related modulation on beta activity of single-trial magnetoencephalography measuring using independent component analysis (ICA), in: Proceedings, 1stInternational IEEE EMBS Confer-ence on Neural Engineering, pp. 396 –398.
Yabe, H., Satio, F., Fukushima, Y., 1993. Median method for detecting endogenous event-related brain potentials. Electroencephalogr. Clin. Neurophysiol. 87, 403– 407.
