- Email: [email protected]
Extracting high frequency oscillatory brain signals from magnetoencephalographic recordings
9th IFAC onBerlin, Biological and Medical Systems Aug. 31 - Symposium Sept. 2, 2015. Germany 9th onBerlin, Biological and Medical Systems Aug.IFAC 31 - Symposium Sept. 2, 2015. Germany 9th IFAC Symposium on Biological and Medicalonline Systems Available at www.sciencedirect.com Aug. 31 - Sept. 2, 2015. Berlin, Germany Aug. 31 - Sept. 2, 2015. Berlin, Germany
ScienceDirect 48-20 (2015) 453–458 ExtractingIFAC-PapersOnLine high frequency oscillatory brain signals from Extracting high frequency oscillatory brain signals from Extracting high frequency oscillatory brain signals from magnetoencephalographic recordings Extracting high frequency oscillatory brain signals from magnetoencephalographic recordings magnetoencephalographic recordings magnetoencephalographic recordings Tilmann H. Sander*, Heriberto Zavala-Fernandez**, Martin Burghoff*,
Tilmann H. Sander*, Heriberto Zavala-Fernandez**, Martin Burghoff*, Yoshinori Uchikawa***, Lutz Trahms* Tilmann Heriberto Zavala-Fernandez**, Martin Tilmann H. H. Sander*, Sander*, Heriberto Zavala-Fernandez**, Martin Burghoff*, Burghoff*, Yoshinori Uchikawa***, Lutz Trahms* Yoshinori Uchikawa***, Lutz Trahms* Yoshinori Uchikawa***, Lutz Trahms* *Physikalisch-Technische Bundesanstalt, 10587 Berlin, Germany (Tel: +30-3481 7436; *Physikalisch-Technische Bundesanstalt, 10587 Berlin, Germany (Tel: +30-3481 7436; e-mail: [email protected]). *Physikalisch-Technische Bundesanstalt, 10587 Berlin, Berlin, Germany Germany (Tel: (Tel: +30-3481 +30-3481 7436; 7436; *Physikalisch-Technische Bundesanstalt, 10587 e-mail: [email protected]). **Department Electrical Engineering, Technical University, 10623 Berlin, Germany. e-mail: [email protected]). e-mail: [email protected]). **Department Electrical Engineering, Technical University, 10623 Berlin, Germany. ***Division**Department of Electronic Electrical and Mechanical Engineering, Tokyo Denki University, Saitama 350-0394, Japan. Engineering, Technical University, 10623 Berlin, Berlin, Germany. Engineering, Technical University, 10623 Germany. ***Division**Department of Electronic Electrical and Mechanical Engineering, Tokyo Denki University, Saitama 350-0394, Japan. ***Division of Electronic and Mechanical Engineering, Tokyo Denki University, Saitama ***Division of Electronic and Mechanical Engineering, Tokyo Denki University, Saitama 350-0394, 350-0394, Japan. Japan. Abstract: Oscillatory short signal bursts of cortical origin with a frequency around 600 Hz can be Abstract: Oscillatory short signal bursts ofThese corticalso origin with a frequency evoked around high 600 Hz can be measured using magnetoencephalography. calledwith somatosensory frequency Abstract: Oscillatory short cortical aa frequency around 600 can Abstract: Oscillatory short signal signal bursts bursts of ofThese corticalso origin origin with frequency evoked around high 600 Hz Hz can be be measured using magnetoencephalography. called somatosensory frequency oscillationsusing (SE-HFO) are induced by an electrical stimulation at the wrist. Up toevoked now only averages over measured magnetoencephalography. These so high frequency measured using magnetoencephalography. Thesestimulation so called called atsomatosensory somatosensory evoked high frequency oscillations (SE-HFO) are induced by an electrical the wrist. Up to now only averages over several thousand stimulations yieldbyananinterpretable result. Here the wrist. spectral and temporal properties of oscillations (SE-HFO) are electrical at Up to only over oscillations (SE-HFO) are induced induced byananinterpretable electrical stimulation stimulation at the the wrist. Upand to now now only averages averages over several thousand stimulations yield result. Here thetemporal spectral temporal properties of SE-HFOs are exploited through epoch concatenationresult. followed by decorrelation to study SEseveral thousand stimulations yield an interpretable Here the spectral and temporal properties of several thousand stimulations yield an interpretable result. Hereby thetemporal spectral decorrelation and temporal to properties of SE-HFOs aretrial exploited through epoch concatenation followed study SEHFO singleare properties. The algorithm is a type offollowed blind source separation and extractstoanstudy SE-HFO SE-HFOs exploited through epoch by decorrelation SESE-HFOs aretrial exploited through epoch concatenation concatenation followed by temporal temporal decorrelation toanstudy SEHFO single properties. The algorithm is a type of blind source separation and extracts SE-HFO component, shows a preferred phase is after sorting the single using aand wave train atan625 Hz as HFO single which trial properties. properties. The algorithm algorithm type of blind blind sourcetrials separation extracts SE-HFO HFO single trial The aa type of source separation extracts SE-HFO component, which showsphase a preferred phase is after trials usingwhich aand wave train atanhave 625 Hz as template. The preferred is not visible aftersorting sortingthe thesingle raw data trials, certainly more component, which shows a preferred phase after sorting the single trials using a wave train at 625 Hz as component, which showsphase a preferred phase after sorting the single trialstrials, usingwhich a wave train at have 625 Hz as template. The preferred is not visible after sorting the raw data certainly more noise. This indicates that the epoch concatenation temporal decorrelation iswhich a powerful toolhave to study template. The preferred phase is not visible after sorting the raw data trials, certainly more template. The preferred phase is not visible after sorting the raw data trials,iswhich certainly more noise. This indicates that the epoch concatenation temporal decorrelation a powerful toolhave to study transient oscillatory signals in multichannel recordings. noise. This indicates that the epoch concatenation temporal decorrelation is a powerful tool to study noise. This indicates that the epoch concatenation temporal decorrelation is a powerful tool to study transient oscillatory signals in multichannel recordings. transient oscillatory signals multichannel recordings. © 2015, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. transient oscillatory signals in in multichannel recordings. Keywords: Magnetoencephalography, independent component analysis, temporal decorrelation, Keywords: Magnetoencephalography, independent component analysis, temporal decorrelation, somatosensory evoked fields, high frequency oscillations, single trials.analysis, temporal decorrelation, Keywords: independent component Keywords: Magnetoencephalography, Magnetoencephalography, independent component somatosensory evoked fields, high frequency oscillations, single trials.analysis, temporal decorrelation, somatosensory evoked fields, high frequency oscillations, somatosensory evoked fields, high frequency oscillations, single single trials. trials. function (PDF) results as the PDF is identically zero between 1. INTRODUCTION function (PDF)wavelets. results asEstimates the PDF isofidentically two SE-HFO PDFs arezero the between basis of function (PDF) results the zero between 1. INTRODUCTION function (PDF)wavelets. results as asEstimates the PDF PDF is isofidentically identically zero between two SE-HFO PDFs are thehas basis of 1. INTRODUCTION independent component analysis of (ICA), which been 1. INTRODUCTION SE-HFO wavelets. Estimates PDFs are the basis of To date somatosensory evoked high frequency oscillations two two SE-HFO component wavelets. Estimates of PDFswhich are thehas basis of independent analysis (ICA), been To date somatosensory evoked high frequency oscillations widely applied to multi-channel data such as MEG and EEG independent component analysis (ICA), which has been (SE-HFO) (Curio et al. (1994); Hashimoto et al. (1996); To date somatosensory evoked high frequency oscillations independent component analysis (ICA), which has been widely applied to multi-channel data such as MEG and EEG To date somatosensory high frequency oscillations (SE-HFO) (Curio al. evoked (1994); Hashimoto (1996); recordings and to functional magnetic resonance imaging (Lee widely applied multi-channel data such and Curio (2004)) are et among the weakest signalset ofal.neuronal (SE-HFO) (Curio et al. Hashimoto et (1996); applied multi-channel dataresonance such as as MEG MEG and EEG EEG recordings and to functional magnetic imaging (Lee (SE-HFO) (Curio et al. (1994); (1994); Hashimoto et ofal. al.neuronal (1996); widely Curio (2004)) are among the weakest signals (1998); Haykin (2000); Hyvarinen et al. (2001); Jung et al. origin, (2004)) that can are be measured using magnetoencephalography recordings and functional magnetic resonance imaging (Lee Curio among the weakest signals of neuronal recordings and functional magnetic resonance imaging (Lee (1998); Haykin (2000); Hyvarinen et al. (2001); Jung et al. Curio (2004)) are among theusing weakest signals of neuronal (2001a); origin, that can be measured magnetoencephalography Roberts et al. (2001); Stone (2002a); Cichocki et al. Haykin (2000); Hyvarinen et al. (2001); Jung (MEG).that These SE-HFOs consist ofmagnetoencephalography a short wavelet of 5-10 (1998); et al. origin, can be measured using (1998); Haykin (2000); Hyvarinen et al. (2001); Jung et (2001a); Roberts et al. (2001); Stone (2002a); Cichocki al. origin, that can be measured using magnetoencephalography (MEG). These SE-HFOs consist of a short wavelet (2002); Makeig etetal. (2004); Choi et(2002a); al. (2005); James et Roberts al. (2001); Stone Cichocki al. ms duration with a maximum amplitude in the range of of 5-10 tens (2001a); et (MEG). These SE-HFOs consist of a short wavelet of 5-10 Robertsetetal. al.(2004); (2001);Choi Stoneet(2002a); Cichocki (2002); al. (2005); James et al. (MEG). These SE-HFOs consist of a short wavelet of ms duration with a maximum amplitude in on the range of 5-10 tens (2001a); (2005)).Makeig (2002); of fT and a frequency around 600 Hz riding top of the well ms duration with amplitude in the range of tens (2002); Makeig Makeig et et al. al. (2004); (2004); Choi Choi et et al. al. (2005); (2005); James James et et al. al. ms duration with aa maximum maximum amplitude in on thetop range of well tens (2005)). of fT and a frequency around 600 Hz riding of the (2005)). known N20m response, which has a latency of 20 ms. The of fT aa frequency around 600 Hz on top of the (2005)). In the literature numerous ICA algorithms have been of fT and and frequency around 600has Hzariding riding on of top20 of ms. the well well known N20m response, which latency The N20m response and the SE-HFOs area induced by an electrical In the literature numerous ICA nature algorithms been known N20m which latency of 20 ms. proposed. Fortunately the sparse of thehave SE-HFOs known N20m response, response, which has has a induced latency by of an 20 electrical ms. The The In the numerous ICA algorithms have been N20m response andwrist the SE-HFOs are stimulation at the (medianus nerve). The SE-HFOs are In the literature literature numerous ICA nature algorithms have been proposed. Fortunately the sparse of the SE-HFOs N20m response and the SE-HFOs are induced by an electrical implies a partitioning ofthethesparse data into a noise andSE-HFOs a signal N20m response andwrist the SE-HFOs are induced bySE-HFOs an electrical proposed. Fortunately nature of the stimulation at the (medianus nerve). The are seen best after averaging several thousand evoked epochs and proposed. Fortunately ofthethesparse nature of the SE-HFOs implies a partitioning data into a noise and a signal stimulation at the wrist (medianus nerve). The SE-HFOs are partitioningofisthe suggested for the application of stimulation at the wrist (medianus nerve). evoked The SE-HFOs are section. implies aaApartitioning data aa noise and seen best after averaging thousand epochs The and suppressing the N20m several response using evoked a bandpass. implies ofisthe data into intofor noise and aa signal signal section. Apartitioning partitioning suggested the application of seen best after averaging several thousand epochs and algorithms such as Infomax (Jung et al. (2001a)) to stimulus seen best afterthe averaging several thousand evoked epochs The and section. A partitioning is suggested for the application of suppressing N20m using a bandpass. extensive literature on the response morphology of averaged SE-HFOs section. A such partitioning is suggested for(2001a)) the application of algorithms as Infomax (Jung et al. to stimulus suppressing the N20m response using a bandpass. The signals. Another partitioning approach is the use of an suppressing the N20m response using a bandpass. The evoked such as (Jung al. to stimulus extensive literature on or theelectroencephalography morphology of averaged SE-HFOs investigated by MEG (EEG) can algorithms algorithms suchAnother as Infomax Infomax (Jung et etapproach al. (2001a)) (2001a)) to use stimulus evoked signals. partitioning is the of an extensive literature on the morphology of averaged SE-HFOs averagedsignals. event related field. In an EEG study is(Porcaro et al. extensive literature on or theelectroencephalography morphology of averaged SE-HFOs Another partitioning approach the an investigated by MEG can evoked be found in Curio (2004). To assess fluctuations of(EEG) SE-HFOs evoked signals. Anotherfield. partitioning approach the use use of of an averaged event related In anofEEG study is(Porcaro et al. investigated by or (EEG) can (2009)) the sequential activation sub-cortical and cortical investigated by MEG MEG or electroencephalography electroencephalography (EEG) can averaged event related field. In an EEG study (Porcaro et al. be found in Curio (2004). To assess fluctuations of SE-HFOs on a single trial (single evoked response) level is of notSE-HFOs possible averaged event related activation field. In anofEEG study (Porcaro et al. (2009)) the sequential sub-cortical and cortical be found in Curio (2004). To assess fluctuations SE-HFOs is demonstrated using sub-cortical the N20 response as be found in Curio (2004). To assess fluctuationsis of SE-HFOs the activation and on a to single (single evoked response) not possible due theirtrial weakness compared to otherlevel signals present in (2009)) (2009)) the issequential sequential activation of of sub-cortical and cortical cortical SE-HFOs the N20 response as on aa to single trial (single evoked level is possible constraint for andemonstrated ICA solution.using on single trial (single compared evoked response) response) level is not not possible SE-HFOs is demonstrated using the N20 response as due their weakness to other signals present in the MEG data. In thiscompared work a to statistical method for the SE-HFOs is andemonstrated constraint for ICA solution.using the N20 response as due to their weakness other signals present in due to theirdata. weakness compared to other signals present in constraint for an ICA solution. the MEG In this work a statistical method for the extraction of single trial SE-HFOs is successfully applied and constraint for anconcatenation ICA solution. temporal decorrelation (ecTD) the data. In this work method for the the epoch the MEG MEG of data. In trial this SE-HFOs work aa statistical statistical method for and the Here extraction single is successfully applied Here the epoch concatenation temporal decorrelation properties of the single trials are derived. extraction single trial SE-HFOs is successfully applied and ICA method (Zavala-Fernandez et al. (2012)) is used(ecTD) on the of the the epoch concatenation temporal decorrelation (ecTD) extraction single trialtrials SE-HFOs is successfully applied and Here properties of single are derived. Here the epoch concatenation temporal decorrelation (ecTD) ICA method (Zavala-Fernandez et al. (2012)) is used on the properties of the single trials are derived. SE-HFO signal section of the raw data. It consists of (Zavala-Fernandez et al. is on properties of the single trials are derived. The isolation of SE-HFO single trials is complicated by their ICA ICA method method (Zavala-Fernandez et raw al. (2012)) (2012)) is used used on the the SE-HFO signal section of the data. It consists of bandpass filtering the raw data between 450 and 1000 Hz, signal of the raw It of The single is complicated by their SE-HFO shortisolation durationof inSE-HFO comparison to trials a typical electrical stimulus SE-HFO signal section section ofdata thebetween raw data. data. Itandconsists consists of bandpass filtering the raw 450 1000 Hz, The isolation of SE-HFO single trials is complicated by their followed filtering by the concatenation of equal450 length windows The isolation ofinSE-HFO single trials is complicated by their bandpass the raw data between and 1000 Hz, short duration comparison to a typical electrical stimulus repetition rate in of comparison 3-9 Hz. ThetoSE-HFOs are easily stimulus seen in bandpass filtering the raw data between 450 and 1000 Hz, followed by the concatenation of equal length windows short duration a typical electrical covering the SE-HFO section of the data. The resulting time short duration in comparison toSE-HFOs a typical electrical stimulus by the of equal length windows repetition are easilyfor seen in followed only 5 msrate of of the3-9 111Hz. msThe interstimulus interval 9 Hz followed by SE-HFO the concatenation concatenation of data. equalThe length windows covering the of the resulting time repetition rate of 3-9 Hz. The SE-HFOs easily in series is subjected to section time-delayed (temporal) decorrelation repetition rate of 3-9 Hz. The SE-HFOs are are easilyforseen seen in covering the SE-HFO section of the data. The resulting time only 5 ms of the 111 ms interstimulus interval 9 Hz stimulation. This means that SE-HFOs are for sparsely covering the SE-HFO section of the (temporal) data. The resulting time series is subjected to time-delayed decorrelation only 5 ms of the 111 ms interstimulus interval 9 Hz (TDD)is(Belouchrani ettime-delayed al. (1997); (temporal) Ziehe et al.decorrelation (1998); de only 5 ms ofThis the 111 ms interstimulus interval for 9 Hz series subjected to stimulation. means that SE-HFOs are sparsely distributed in measured raw data corresponding to nonseries is subjected to time-delayed (temporal) decorrelation (TDD) (Belouchrani et al.TDD (1997); et al. (1998); stimulation. This SE-HFOs are Cheveigne et al. (2014)). was Ziehe successfully applied de in stimulation. This means meansrawthat that are sparsely sparsely (Belouchrani et (1997); Ziehe et de distributed dataSE-HFOs corresponding to nonstationarity: in To measured describe statistically the SE-HFO signal on a (TDD) (TDD) (Belouchrani et al. al.TDD (1997); Ziehe et al. al. (1998); (1998); Cheveigne et al. (2014)). was successfully applied de in distributed in measured raw data corresponding to nonthe study ofetN20m responses (Tang etsuccessfully al. (2002); Kishida et distributed in measured raw data the corresponding to nonstationarity: To describe statistically SE-HFO signal on a Cheveigne al. (2014)). TDD was applied 10 ms timeToscale a time dependent probability density Cheveigne etN20m al. (2014)). TDD was etsuccessfully applied in in the study of responses (Tang al. (2002); Kishida et stationarity: describe statistically the SE-HFO signal on a (2003); Tang et al. (2005); Sander et al. al. (2005)). stationarity: describe statistically the SE-HFO signal on a al. the study of N20m responses (Tang et (2002); Kishida 10 ms timeToscale a time dependent probability density the study ofTang N20m responses (Tang et (2002); Kishida et et al. (2003); et al. (2005); Sander et al. (2005)). 10 ms time scale a time dependent probability density 10 ms time scale a time dependent probability density al. (2003); Tang et al. (2005); Sander et al. (2005)). al. (2003); Tang et al. (2005); Sander et al. (2005)).
Copyright 2015 453 Hosting by Elsevier Ltd. All rights reserved. 2405-8963© 2015,IFAC IFAC (International Federation of Automatic Control) Copyright ©©2015 IFAC 453 Copyright © 2015 453Control. Peer review underIFAC responsibility of International Federation of Automatic Copyright © 2015 IFAC 453 10.1016/j.ifacol.2015.10.182
9th IFAC BMS 454 Tilmann H. Sander et al. / IFAC-PapersOnLine 48-20 (2015) 453–458 Aug. 31 - Sept. 2, 2015. Berlin, Germany
2. METHODS 2.1 Concatenated data and temporal decorrelation Given a superposition of n sources of the form ai si(t), where ai is a time independent field map of the m channel MEG system and si(t) describes the time dependence of the source, the measured data vector x(t) can be decomposed by a demixing matrix W: W x(t) = u(t).
(1)
Here the ui(t) are identical to the si(t) if disregarding an indeterminacy in scale and ordering. To estimate a W blindly ICA is frequently applied and n=m is assumed. One ICA method suitable for MEG data is the time-delayed decorrelation (TDD, TDSEP, SOBI) algorithm (Belouchrani et al. (1997); Ziehe et al. (1998); Congedo et al. (2008)), which is often classified as a blind source separation type algorithm as it relies solely on second-order statistics. In TDD a set of time-delayed covariance matrices is calculated, these matrices are symmetrised, and subsequently approximately diagonalised to obtain an estimate of W. As the SE-HFOs have a limited duration and a specific frequency the ecTD algorithm (Zavala-Fernandez et al. (2012)) derived from TDD seems suitable to extract them. Firstly the data are bandpass filtered between 450 and 1000 Hz. Then equal length windows covering the duration of the SE-HFO wavelet are concatenated to form a new continuous time series. Finally time-delayed covariance matrices of the concatenated time series are calculated and processed as in standard TDD.
Fig. 1. Illustration of input data for the temporal decorrelation ICA. a) TDD: Two time series x and y from a multi-channel measurement are used in the time-delayed covariance calculation C(x,y,τmax). b) Epoch concatenation TD: Concatenation of equal length time windows wi of x and y yields new time series xc and yc used for C(xc,yc,τmax).
To contrast ecTD with TDD the input data for the maximum delay τmax covariance matrix are illustrated in Fig. 1 for both algorithms. In Fig. 1a) two time series x and y of a multichannel measurement are sketched. The data range 0...T of x is correlated with the data range τmax...T+τmax of y (ranges indicated by red bar) to obtain the matrix element C(x,y,τmax) of TDD.
noise level. This avoids steps in the concatenated time series. For ecTD the τmax covariance matrix C(xc,yc,τmax) is calculated from concatenated time series such as xc and yc for the data ranges indicated in red. The SE-HFOs shown in Fig. 1b) are exaggerated and cannot be observed in real data.
Due to the wavelet like nature of SE-HFOs they can be modeled as s(t) = A(t) sin(2 π fHF t + φ)
(2)
The C(xc,yc,τk) matrix of the ecTD algorithm for M windows of length wlen is given by
with envelope A(t), frequency fHF, and phase φ. In this model f is about 600 Hz and the envelope A(t) has a support of less than 10 ms centered on the N20m peak. This model is motivated by physiological knowledge derived from electrical recordings, which show that the SE-HFOs are always associated with the N20m peak in a systematic fashion (Klostermann et al. (1999)).
Cij(xc,yc,τk) = Σ{l} xi(tl) xj(tl+ τk),
(3)
where the set of indices is {l} = w1,start,...,w1,start + wlen, w2,start,...,w2,start + wlen,wM,start,...,wM,start + wlen and contains M* wlen elements. The ecTD algorithm requires choosing three parameters: The window start wi,start relative to the stimulus, the window length wlen, and the maximum delay τmax. The window start and the window length are fairly obvious in the case of SEHFOs as the window should contain the complete SE-HFO wavelet. The length of the SE-HFO can be estimated from a conventional average. The choice of τmax is largely empirical,
The raw SE-HFO input data for ecTD are sketched in Fig. 1b) as a new concatenated time series xc is created from x using the data in the windows wi of length wlen indicated in blue. Any discontinuity between the windows will be obscured in the noise as the length of the windows extends beyond the SE-HFOs, i.e., the range where A(t) is above 454
9th IFAC BMS Aug. 31 - Sept. 2, 2015. Berlin, Germany Tilmann H. Sander et al. / IFAC-PapersOnLine 48-20 (2015) 453–458
455
an accepted rule of thumb is to set a desirable frequency resolution f and set τmax= 1/f.
but it was found that ecTD works best for 1/τmax >> wlen. Using a time-delay larger than the window length implies similar data and noise content in each window, i.e., stationarity. A rule of thumb for TDD states that the frequency resolution is 1/τmax (Sander et al. (2005)).
Once single trial SE-HFOs have been extracted they can be characterised using (2). Set values for fHF and a time range tstart ... tstop are chosen and then time independent average amplitudes Ai and phases φi are calculated for each trial (Delorme et al. (2003)). Only the three central oscillations of the wavelet are used to estimate Ai. Subsequently the single trials can be sorted by amplitude or phase or by a sequence of both. For visualization after sorting a single trial (ST) plot (Jung et al. (2001a); Tang et al. (2002)) will be used. In the ST plot multiple sections of a time series are displayed as horizontal scan lines stacked vertically with the function value coded as a grey shade between black and white.
Fig. 2. Conventionally averaged MEG result for an electrical stimulation at the left medianus nerve. The map sequence shows at 21 ms a dipolar pattern known as the N20m. Channel CH44 has the highest peak amplitude and small variations in the shape of the peak are the SE-HFOs directly visible in Fig. 3 after bandpass filtering. 2.2 Measurements and data processing In one healthy volunteer MEG was recorded using a 125channel whole head MEG system made by Eagle Technology (Yokogawa MEG Vision (1998)). To elicit a SEF N20m an electrical pulse of 100 μs duration and 9 Hz repetition rate was applied to the left medianus nerve with a strength of 12 mA selected to be above the motor threshold. The measurement lasted 670 s and the data were sampled at 5 kHz. All measurements were performed in accordance with the declaration of Helsinki. Data from other subjects showed similar results.
Fig. 3. Average over the highpass filtered data shown as a sequence of maps together with the result for a single channel (top part) and the HF-ecTD component map with the associated averaged time series (bottom part). The oscillations in both time series show the same phase and envelope and the SE-HFOs have a frequency of f ≈ 625 Hz.
Before the ecTD calculation the multi-channel raw data set was digitally filtered using a 450 to 1000 Hz bandpass filter with -10 dB damping points at 400 and 1100 Hz and no phase shifts in the passband. After digital filtering the epochs were concatenated to create a new data set (see section 2.1) using wstart = 16 ms post stimulus and wlen = 13 ms. This choice of wstart and wlen is valid for any other dataset of this type as the N20m occurs subject independent at 20 ms with a variation of at most 2 ms and therefore this choice covers the SE-HFOs riding on top of the N20m.
3. RESULTS 3.1 Extraction of SE-HFO signal The conventional N20m result is given in Fig. 2, which shows a map sequence from the average over 6000 epochs of recorded data and the time series of channel CH44, which has the highest peak amplitude at 20 ms. A 5 ms baseline starting at 5 ms after the stimulus was subtracted.
The set of time delayed covariance matrices was chosen empirically to be {τi} = {2, 4, 6,..., 2000 ms}, which requires to diagonalise 1000 matrices. The largest τi corresponds to a frequency of f = 1/τmax = 1/2000 ms = 0.5 Hz. A theoretical foundation for the choice of τmax for TDD is not known and
The map at t = 0 s is the stimulus artifact due to the electrical current pulse, at 16 ms no activity is seen, and a clear dipolar pattern is observed in the map at 21 ms. Later activity at 24 ms is more complex. The dipolar pattern at 21 ms is on the 455
9th IFAC BMS 456 Tilmann H. Sander et al. / IFAC-PapersOnLine 48-20 (2015) 453–458 Aug. 31 - Sept. 2, 2015. Berlin, Germany
right side of the map contralateral to the left stimulated side (orientation of map: top of map points to nose, left and right of map correspond to left and right hemisphere). A dipolar pattern as shown in the map at 21 ms was successfully approximated by an ECD source in the cortex for the same whole head MEG system (Sander et al. (2007)). To extract the SE-HFOs the raw data have to be filtered using the 450-1000 Hz bandpass mentioned previously (see section 2.2) and then averaged. The SE-HFOs appear then as a short wavelet around 21 ms in the average. This is shown in the top part of Fig. 3, where a map sequence and the single channel of the averaged filtered data are given. The SE-HFOs can be seen in Fig. 2 as small peak shape variations in the same channel (CH44) before filtering. The maps in Fig. 3 appear much noisier compared to the maps in Fig. 2, which is not surprising noting the different contour step for the maps in the Figs. 2 and 3 and the well known small amplitude of the SE-HFOs. The maps in Fig. 3 show a dipolar structure with a clear resemblance to the N20m map in Fig. 2. The time points of the maps were chosen to coincide with the oscillation extrema in the time series shown below the maps. This time series shows the channel Ch 44, which had the strongest oscillation amplitude. The location of this channel is indicated in the map. As the maps are aligned with the oscillation extrema the sign of the map reverses with a periodicity of 1.6 ms (= 625 Hz), the periodicity indicated for the time series. The first peak of the SE-HFOs appears at 17 ms, which is consistent with the literature (Curio et al. (1994); Hashimoto et al. (1996)).
Fig. 4. Comparison of power spectra of unaveraged data. Top: Bandpass filtered data of Ch 44. Middle: Concatenated bandpass filtered data of of Ch 44 (wlen = 13 ms). Bottom: The time series of the HF-ecTD component from Fig. 3. The peaks marked by red circles in the middle panel dominate the spectrum of the HF-ecTD component.
the concatenated data. The occurrence of several peaks in the spectrum of HF-ecTD can be explained as an amplitude modulation process. Assuming a single frequency f modulated by an envelope fm peaks at f+fm and f-fm result. The concatenated SE-HFOs can be viewed as such a type of signal with f ≈ 625 Hz and fM = 77 Hz. As the SE-HFOs do not have a simple sinusoidal envelope more peaks result in the spectrum. It was verified that changing the concatenation window length changes the peak separation. The ecTD algorithm apparently succeeded in isolating the SE-HFO signal due to its periodic signature above noise level in the concatenated signal. In agreement with that no component related to the SE-HFOs was obtained by applying TDD to the full length band pass filtered data using similar time-delay parameters as used in the ecTD calculation.
As an alternative method of extracting the SE-HFOs the ecTD algorithm was applied to the bandpass filtered data using the parameters described in the previous section. Only one of the 125 resulting component maps had a dipolar pattern, which is a sign of a brain source (Delorme et al. (2012)). This component map is shown in the bottom part of Fig. 3 and labeled as HF-ecTD. The average over its time series is similar to the Ch 44 time series as can be seen in the middle of Fig. 3 and therefore the component HF-ecTD is attributed to the SE-HFOs. All other ecTD components (not shown) have an irregular map typical of noise components. The HF-ecTD average time series covers only the range from 15 to 28 ms as this was the concatenation window for the ecTD algorithm. The HF-ecTD component map is similar to the map at 23 ms from the simple average (Fig. 3, top). The similarity of the two maps can be quantified by an angle α defined as cos α = a*b / ( |a|*|b| ) and α = 26o results, which means a high degree of similarity given the additional noise visible in the maps.
3.2 Single trial analysis The SE-HFOs were analysed subsequently using the model in (2). The parameters Ai and phase φi of each trial were estimated for fHF = 625 Hz and tstart ... tstop covering a 4.8 ms window around t = 23.6 ms post-stimulus. The choice of tstart and tstop was motivated by the envelope of CH44 in Fig. 2, which shows three peaks at the same amplitude followed by two peaks with a higher amplitude. The existence of two regimes in the SE-HFOs is well known (Klostermann et al. (1999); Haueisen et al. (2007)). This choice of tstart and tstop might be not suitable for other datasets and similarly the averaged SE-HFOs would need to be inspected.
An important property of the HF-ecTD component is presented in Fig. 4, which compares the power spectra of the unaveraged data for different processing stages. The spectrum of the bandpass filtered Ch 44 data does not have specific features, whereas four peaks separated by f = 1.0/ wlen = 77 Hz (= 1.0/13 ms) are observed above the background in the spectrum of the concatenated data (peaks indicated by red circles). Very distinct peaks are visible in the spectrum of HF-ecTD at the same positions as the peaks in
Then the trials were sorted by Ai (equation (2)), split into four groups, and sorted again within each group by phase φi. The result is shown in Fig. 5, where the ST plots (see section 2.2) of the high and low amplitude groups (top, middle) are shown together with the subaverages of the four groups 456
9th IFAC BMS Aug. 31 - Sept. 2, 2015. Berlin, Germany Tilmann H. Sander et al. / IFAC-PapersOnLine 48-20 (2015) 453–458
(bottom; orange = low amplitude, dark blue = high amplitude) for the time series of Ch 44 and of component HF-ecTD. For display purposes non-overlapping groups of ten trials were averaged and an 11 ms window is shown. The sign is changed deliberately in Fig. 5 between Ch 44 and HFecTD to facilitate a comparison of the ST plots.
457
Nevertheless the total average in Fig. 3 is very similar for the two ensembles of trials, i.e., both ensembles contain the same amount of signal power. The observation of a preferred phase in the HF-ecTD-ST plot is only possible through the separation of noise and signal in the ecTD calculation. Note that the four subgroup averages in Fig. 5 (bottom) of Ch 44 appear to have a higher noise level compared to those of HFecTD in agreement with a dominant noise removal effect. Only a single dipolar pattern, component HF-ecTD, was observed in the ecTD result despite the possibility to detect non-orthogonal sources using an ICA-type method. The oscillatory signal in the four SE-HFOs subaverages in Fig. 5 (bottom) is almost identical outside the sorting window, i.e., for t < 21 ms. Apparently the double sorting is relevant only for the 5 ms sorting window. This indicates that the weaker oscillations between 16 and 21 ms could originate from a different process compared to the stronger oscillations between 21 and 26 ms. This might be related to the early and late (before/after N20m) SE-HFO regimes observable in EEG and weaker in MEG data (Klostermann et al. (1999); Haueisen et al. (2007)). 4. CONCLUSION The successful extracting of an ensemble of SE-HFOs single trials by ecTD shows that ecTD is useful to extract transient oscillatory signals in multichannel recordings such as MEG. The ensemble extracted by ecTD is much better suited for further processing compared to the original ensemble as can be seen from the single trial analysis. Although this methodology was demonstrated on a single dataset it will be applicable to study group effects, because the SE-HFOs associated with the N20m can be observed in most subjects. The application of standard TDD (SOBI, TDSEP) to N20m data (Tang et al. (2002); Kishida et al. (2003); Tang et al. (2005); Sander et al. (2005)) is directly related to the present situation. The signal repetition and partitioning created through concatenation in ecTD is implicit in raw N20m data as the stimulation is periodic and the N20m neuronal response is strong enough enough to yield a distinct spectral pattern. For the N20m the temporal support has a sufficient length for temporal decorrelation to work, but for SE-HFOs the temporal extension is too small and concatenation is needed.
Fig. 5. Sorted ST plots for channel Ch 44 (left) and the time series of the HF-ecTD component (right). The trials were sorted by amplitude Ai and phase φi in four groups. Of the four groups the ST plots of the highest and lowest amplitude groups are shown (top, middle) together with the individual subaverages for the four groups (bottom).
Comparing upper and lower ST plot in Fig. 5 it is evident that the single trial SE-HFOs cover a wide range of amplitudes. Comparing the high amplitude ST plots of Ch 44 and HFecTD shows that a preferential phase exists in the single trials of HF-ecTD as most of the single trials have a peak aligned with the black vertical line added as a guide to the eye. The peaks of the HF-ecTD single trials form a sigmoid line, it is almost straight for Ch 44 corresponding to a uniformly distributed phase. The average over an oscillation with a uniformly distributed phase is zero for an infinite number of trials and mainly the amplitude distribution of the trials in Ch 44 leads to a non-zero average.
Clearly the known modulations of the SE-HFOs (Curio (2004)) should be investigated by extracting single trials using ecTD on a group of subjects, e.g., is the decay of the SE-HFO amplitude during sleep due to phase randomization or due to a decrease of amplitude? Newer studies (Fedele et al. (2015); Leiken et al. (2014)) have measured noninvasively electro-physiological brain signals around 1 kHz and ecTD might be suitable to get more insights into these data as well. ACKNOWLEGDMENTS This work is supported by BMBF grant 01 GO 0208 BNIC (Berlin Neuroimaging Centre) and DAAD scholarship A0421558.
Related to this are the grey scale bars on top of the Ch 44and HF-ecTD-ST plots showing that the single trials in Ch 44 have typically twice the amplitude compared to HF-ecTD. 457
9th IFAC BMS 458 Tilmann H. Sander et al. / IFAC-PapersOnLine 48-20 (2015) 453–458 Aug. 31 - Sept. 2, 2015. Berlin, Germany
REFRENCES
Kishida K., Kato K., Shinosaki K., Ukai S. (2003). Blind identification of SEF dynamics from MEG data by using decorrelation method of ICA. Proc. of the 4th Intl. Symposium on ICA nd BSS, Nara, 185-190. Klostermann F., Nolte G., Curio G. (1999). Multiple generators of 600 Hz wavelets in human SEP unmasked by varying stimulus rates. Neuroreport, 10, 1625-1629. Lee T.W. (1998). Independent Component Analysis - Theory and Applications. Kluwer Academic Publishers, Boston. Leiken K., Xiang J., Zhang F., Shi J., Tang L., Liu H., Wang X. (2014). Magnetoencephalography detection of highfrequency oscillations in the developing brain. Front. Hum. Neurosci., 8, 969. Makeig S., Debener S., Onton J., Delorme A. (2002). Mining event-related brain dynamics. TRENDS in Cogn. Sciences, 8, 204-210. Porcaro C., Coppola G., DiLorenzo G., Zappasodi F., Siracusano A., Pierelli F., Rossini P.M., Tecchio F., Seri S. (2009). Hand somatosensory subcortical and cortical sources assessed by functional source separation: An EEG study. Hum. Brain Mapp., 30, 660-674. Roberts S., Everson R. (eds.) (2001). Independent Component Analysis - Principles and Practice. Cambridge University Press, Cambridge. Sander T.H., Burghoff M., Curio G., Trahms L. (2005). Single evoked somatosensory MEG responses extracted by time delayed decorrelation. IEEE Trans. on Signal Processing, 53, 3384-3392. Sander T.H, Burghoff M., Van Leeuwen P., Trahms L. (2007) Application of decorrelation independent component analysis to biomagnetic multi-channel measurements. Biomedizinische Technik 52: 130–136. Stone J.V. (2002). Independent component analysis: An introduction. TRENDS in Cogn. Sciences, 6, 59-64. Tang A.C., Pearlmutter B.A., Malaszenko N.A., Phung D.B. (2002). Independent components of magnetoencephalography: Single-trial response onset times. Neuroimage, 17, 1773-1789. Tang A.C., Sutherland M.T., McKinney C.J. (2005). Validation of SOBI components from high-density EEG. Neuroimage, 25, 539-553. Yokogawa MEG Vision (1998). Distributed by YOKOGAWA / Japan. http://www.yokogawa.com/us/technical-library/whitepapers/megvision-magnetoencephalograph-system-andits-applications.htm (accessed 23.1.2015). Zavala-Fernandez H., Orglmeister R., Trahms L., Sander T.H. (2012). Identification enhancement of auditory evoked potentials in EEG by epoch concatenation and temporal decorrelation. Comput Methods Programs Biomed. 108, 1097-105. Ziehe A., Müller K.-R. (1998). TDSEP - an efficient algorithm for blind separation using time structure. In L. Niklasson et al. (eds.), Proceedings of the 8th ICANN, 675-680. Springer.
Belouchrani A., Abed-Meraim K., Cardoso J.-F., Moulines E. (1997). A blind source separation technique based on second-order statistics. IEEE Trans. Sig. Proc., 45, 434444. Choi S., Cichocki A., Park H.-M., Lee S.-Y. (2005). Blind source separation and independent component analysis: A review. Neural Inf. Proc., 6, 1-57. Cichocki A., S.Amari S. (2002). Adaptive Blind Signal and Image Processing . John Wiley and Sons, New York. Congedo M., Gouy-Pailler C., Jutten C. (2008). On the blind source separation of human electroencephalogram by approximate joint diagonalization of second order statistics. Clin. Neurophysiol., 119, 2677-2686. Curio G., Mackert B.-M., Burghoff M., Koetitz R., AbrahamFuchs K., Härer W. (1994). Localization of evoked neuromagnetic 600Hz activity in the cerebral somatosensory system, Elec. Clin. Neurophysiol., 91, 483-487. Curio G. (2004). Ultrafast EEG Activities. In E. Niedermeyer and L. da Silva (eds.), Electroencephalography - Basic Principles, Clinical Applications, and Related Fields, 495-504. Lippincott, Williams and Wilkins, Phialdelphia. de Cheveigné A., Parra L.C. (2014). Joint decorrelation, a versatile tool for multichannel data analysis. Neuroimage, 98, 487-505. Delorme A., Makeig S. (2003). EEGLAB: An open source toolbox for analysis of single-trial EEG dynamics including independent component analysis. J. of Neurosc. Meth., 134, 9-21. Delorme A., Palmer J., Onton J., Oostenveld R., Makeig S. (2012). Independent EEG sources are dipolar. PLoS One, 7, e30135. Fedele T, Scheer HJ, Burghoff M, Curio G, Körber R. (2015). Ultra-low-noise EEG/MEG systems enable bimodal non-invasive detection of spike-like human somatosensory evoked responses at 1 kHz. Physiol. Meas., 36, 357-368. Hashimoto I., Mashiko T., Imada T. (1996). Somatic evoked high-frequency magnetic oscillations reflect activity of inhibitory interneurons in the human somatosensory cortex. Elec. Clin. Neurophysiol., 100, 189-203. Haueisen J., Leistritz L., Süsse T., Curio G., H.Witte H. (2007). Identifying mutual information transfer in the brain with differential-algebraic modeling: Evidence for fast oscillatory coupling between cortical somatosensory areas 3b and 1. Neuroimage, 37, 130-136. Haykin S. (ed.) (2000). Unsupervised Adaptive Filtering. Volume I: Blind Source Separation. John Wiley and Sons, New York. Hyvärinen A., Karhunen J., Oja E. (2001). Independent Component Analysis. John Wiley and Sons, New York. James C.J., Hesse C.W. (2005). Independent component analysis for biomedical signals. Physiol. Meas., 26, R15R39. Jung T.P., Makeig S., McKeown M.J., Bell A.J., Lee T.W., Sejnowski T. (2001). Imaging brain dynamics using independent component analysis. Proc. of IEEE, 89, 1107-1122.
458
ScienceDirect 48-20 (2015) 453–458 ExtractingIFAC-PapersOnLine high frequency oscillatory brain signals from Extracting high frequency oscillatory brain signals from Extracting high frequency oscillatory brain signals from magnetoencephalographic recordings Extracting high frequency oscillatory brain signals from magnetoencephalographic recordings magnetoencephalographic recordings magnetoencephalographic recordings Tilmann H. Sander*, Heriberto Zavala-Fernandez**, Martin Burghoff*,
Tilmann H. Sander*, Heriberto Zavala-Fernandez**, Martin Burghoff*, Yoshinori Uchikawa***, Lutz Trahms* Tilmann Heriberto Zavala-Fernandez**, Martin Tilmann H. H. Sander*, Sander*, Heriberto Zavala-Fernandez**, Martin Burghoff*, Burghoff*, Yoshinori Uchikawa***, Lutz Trahms* Yoshinori Uchikawa***, Lutz Trahms* Yoshinori Uchikawa***, Lutz Trahms* *Physikalisch-Technische Bundesanstalt, 10587 Berlin, Germany (Tel: +30-3481 7436; *Physikalisch-Technische Bundesanstalt, 10587 Berlin, Germany (Tel: +30-3481 7436; e-mail: [email protected]). *Physikalisch-Technische Bundesanstalt, 10587 Berlin, Berlin, Germany Germany (Tel: (Tel: +30-3481 +30-3481 7436; 7436; *Physikalisch-Technische Bundesanstalt, 10587 e-mail: [email protected]). **Department Electrical Engineering, Technical University, 10623 Berlin, Germany. e-mail: [email protected]). e-mail: [email protected]). **Department Electrical Engineering, Technical University, 10623 Berlin, Germany. ***Division**Department of Electronic Electrical and Mechanical Engineering, Tokyo Denki University, Saitama 350-0394, Japan. Engineering, Technical University, 10623 Berlin, Berlin, Germany. Engineering, Technical University, 10623 Germany. ***Division**Department of Electronic Electrical and Mechanical Engineering, Tokyo Denki University, Saitama 350-0394, Japan. ***Division of Electronic and Mechanical Engineering, Tokyo Denki University, Saitama ***Division of Electronic and Mechanical Engineering, Tokyo Denki University, Saitama 350-0394, 350-0394, Japan. Japan. Abstract: Oscillatory short signal bursts of cortical origin with a frequency around 600 Hz can be Abstract: Oscillatory short signal bursts ofThese corticalso origin with a frequency evoked around high 600 Hz can be measured using magnetoencephalography. calledwith somatosensory frequency Abstract: Oscillatory short cortical aa frequency around 600 can Abstract: Oscillatory short signal signal bursts bursts of ofThese corticalso origin origin with frequency evoked around high 600 Hz Hz can be be measured using magnetoencephalography. called somatosensory frequency oscillationsusing (SE-HFO) are induced by an electrical stimulation at the wrist. Up toevoked now only averages over measured magnetoencephalography. These so high frequency measured using magnetoencephalography. Thesestimulation so called called atsomatosensory somatosensory evoked high frequency oscillations (SE-HFO) are induced by an electrical the wrist. Up to now only averages over several thousand stimulations yieldbyananinterpretable result. Here the wrist. spectral and temporal properties of oscillations (SE-HFO) are electrical at Up to only over oscillations (SE-HFO) are induced induced byananinterpretable electrical stimulation stimulation at the the wrist. Upand to now now only averages averages over several thousand stimulations yield result. Here thetemporal spectral temporal properties of SE-HFOs are exploited through epoch concatenationresult. followed by decorrelation to study SEseveral thousand stimulations yield an interpretable Here the spectral and temporal properties of several thousand stimulations yield an interpretable result. Hereby thetemporal spectral decorrelation and temporal to properties of SE-HFOs aretrial exploited through epoch concatenation followed study SEHFO singleare properties. The algorithm is a type offollowed blind source separation and extractstoanstudy SE-HFO SE-HFOs exploited through epoch by decorrelation SESE-HFOs aretrial exploited through epoch concatenation concatenation followed by temporal temporal decorrelation toanstudy SEHFO single properties. The algorithm is a type of blind source separation and extracts SE-HFO component, shows a preferred phase is after sorting the single using aand wave train atan625 Hz as HFO single which trial properties. properties. The algorithm algorithm type of blind blind sourcetrials separation extracts SE-HFO HFO single trial The aa type of source separation extracts SE-HFO component, which showsphase a preferred phase is after trials usingwhich aand wave train atanhave 625 Hz as template. The preferred is not visible aftersorting sortingthe thesingle raw data trials, certainly more component, which shows a preferred phase after sorting the single trials using a wave train at 625 Hz as component, which showsphase a preferred phase after sorting the single trialstrials, usingwhich a wave train at have 625 Hz as template. The preferred is not visible after sorting the raw data certainly more noise. This indicates that the epoch concatenation temporal decorrelation iswhich a powerful toolhave to study template. The preferred phase is not visible after sorting the raw data trials, certainly more template. The preferred phase is not visible after sorting the raw data trials,iswhich certainly more noise. This indicates that the epoch concatenation temporal decorrelation a powerful toolhave to study transient oscillatory signals in multichannel recordings. noise. This indicates that the epoch concatenation temporal decorrelation is a powerful tool to study noise. This indicates that the epoch concatenation temporal decorrelation is a powerful tool to study transient oscillatory signals in multichannel recordings. transient oscillatory signals multichannel recordings. © 2015, IFAC (International Federation of Automatic Control) Hosting by Elsevier Ltd. All rights reserved. transient oscillatory signals in in multichannel recordings. Keywords: Magnetoencephalography, independent component analysis, temporal decorrelation, Keywords: Magnetoencephalography, independent component analysis, temporal decorrelation, somatosensory evoked fields, high frequency oscillations, single trials.analysis, temporal decorrelation, Keywords: independent component Keywords: Magnetoencephalography, Magnetoencephalography, independent component somatosensory evoked fields, high frequency oscillations, single trials.analysis, temporal decorrelation, somatosensory evoked fields, high frequency oscillations, somatosensory evoked fields, high frequency oscillations, single single trials. trials. function (PDF) results as the PDF is identically zero between 1. INTRODUCTION function (PDF)wavelets. results asEstimates the PDF isofidentically two SE-HFO PDFs arezero the between basis of function (PDF) results the zero between 1. INTRODUCTION function (PDF)wavelets. results as asEstimates the PDF PDF is isofidentically identically zero between two SE-HFO PDFs are thehas basis of 1. INTRODUCTION independent component analysis of (ICA), which been 1. INTRODUCTION SE-HFO wavelets. Estimates PDFs are the basis of To date somatosensory evoked high frequency oscillations two two SE-HFO component wavelets. Estimates of PDFswhich are thehas basis of independent analysis (ICA), been To date somatosensory evoked high frequency oscillations widely applied to multi-channel data such as MEG and EEG independent component analysis (ICA), which has been (SE-HFO) (Curio et al. (1994); Hashimoto et al. (1996); To date somatosensory evoked high frequency oscillations independent component analysis (ICA), which has been widely applied to multi-channel data such as MEG and EEG To date somatosensory high frequency oscillations (SE-HFO) (Curio al. evoked (1994); Hashimoto (1996); recordings and to functional magnetic resonance imaging (Lee widely applied multi-channel data such and Curio (2004)) are et among the weakest signalset ofal.neuronal (SE-HFO) (Curio et al. Hashimoto et (1996); applied multi-channel dataresonance such as as MEG MEG and EEG EEG recordings and to functional magnetic imaging (Lee (SE-HFO) (Curio et al. (1994); (1994); Hashimoto et ofal. al.neuronal (1996); widely Curio (2004)) are among the weakest signals (1998); Haykin (2000); Hyvarinen et al. (2001); Jung et al. origin, (2004)) that can are be measured using magnetoencephalography recordings and functional magnetic resonance imaging (Lee Curio among the weakest signals of neuronal recordings and functional magnetic resonance imaging (Lee (1998); Haykin (2000); Hyvarinen et al. (2001); Jung et al. Curio (2004)) are among theusing weakest signals of neuronal (2001a); origin, that can be measured magnetoencephalography Roberts et al. (2001); Stone (2002a); Cichocki et al. Haykin (2000); Hyvarinen et al. (2001); Jung (MEG).that These SE-HFOs consist ofmagnetoencephalography a short wavelet of 5-10 (1998); et al. origin, can be measured using (1998); Haykin (2000); Hyvarinen et al. (2001); Jung et (2001a); Roberts et al. (2001); Stone (2002a); Cichocki al. origin, that can be measured using magnetoencephalography (MEG). These SE-HFOs consist of a short wavelet (2002); Makeig etetal. (2004); Choi et(2002a); al. (2005); James et Roberts al. (2001); Stone Cichocki al. ms duration with a maximum amplitude in the range of of 5-10 tens (2001a); et (MEG). These SE-HFOs consist of a short wavelet of 5-10 Robertsetetal. al.(2004); (2001);Choi Stoneet(2002a); Cichocki (2002); al. (2005); James et al. (MEG). These SE-HFOs consist of a short wavelet of ms duration with a maximum amplitude in on the range of 5-10 tens (2001a); (2005)).Makeig (2002); of fT and a frequency around 600 Hz riding top of the well ms duration with amplitude in the range of tens (2002); Makeig Makeig et et al. al. (2004); (2004); Choi Choi et et al. al. (2005); (2005); James James et et al. al. ms duration with aa maximum maximum amplitude in on thetop range of well tens (2005)). of fT and a frequency around 600 Hz riding of the (2005)). known N20m response, which has a latency of 20 ms. The of fT aa frequency around 600 Hz on top of the (2005)). In the literature numerous ICA algorithms have been of fT and and frequency around 600has Hzariding riding on of top20 of ms. the well well known N20m response, which latency The N20m response and the SE-HFOs area induced by an electrical In the literature numerous ICA nature algorithms been known N20m which latency of 20 ms. proposed. Fortunately the sparse of thehave SE-HFOs known N20m response, response, which has has a induced latency by of an 20 electrical ms. The The In the numerous ICA algorithms have been N20m response andwrist the SE-HFOs are stimulation at the (medianus nerve). The SE-HFOs are In the literature literature numerous ICA nature algorithms have been proposed. Fortunately the sparse of the SE-HFOs N20m response and the SE-HFOs are induced by an electrical implies a partitioning ofthethesparse data into a noise andSE-HFOs a signal N20m response andwrist the SE-HFOs are induced bySE-HFOs an electrical proposed. Fortunately nature of the stimulation at the (medianus nerve). The are seen best after averaging several thousand evoked epochs and proposed. Fortunately ofthethesparse nature of the SE-HFOs implies a partitioning data into a noise and a signal stimulation at the wrist (medianus nerve). The SE-HFOs are partitioningofisthe suggested for the application of stimulation at the wrist (medianus nerve). evoked The SE-HFOs are section. implies aaApartitioning data aa noise and seen best after averaging thousand epochs The and suppressing the N20m several response using evoked a bandpass. implies ofisthe data into intofor noise and aa signal signal section. Apartitioning partitioning suggested the application of seen best after averaging several thousand epochs and algorithms such as Infomax (Jung et al. (2001a)) to stimulus seen best afterthe averaging several thousand evoked epochs The and section. A partitioning is suggested for the application of suppressing N20m using a bandpass. extensive literature on the response morphology of averaged SE-HFOs section. A such partitioning is suggested for(2001a)) the application of algorithms as Infomax (Jung et al. to stimulus suppressing the N20m response using a bandpass. The signals. Another partitioning approach is the use of an suppressing the N20m response using a bandpass. The evoked such as (Jung al. to stimulus extensive literature on or theelectroencephalography morphology of averaged SE-HFOs investigated by MEG (EEG) can algorithms algorithms suchAnother as Infomax Infomax (Jung et etapproach al. (2001a)) (2001a)) to use stimulus evoked signals. partitioning is the of an extensive literature on the morphology of averaged SE-HFOs averagedsignals. event related field. In an EEG study is(Porcaro et al. extensive literature on or theelectroencephalography morphology of averaged SE-HFOs Another partitioning approach the an investigated by MEG can evoked be found in Curio (2004). To assess fluctuations of(EEG) SE-HFOs evoked signals. Anotherfield. partitioning approach the use use of of an averaged event related In anofEEG study is(Porcaro et al. investigated by or (EEG) can (2009)) the sequential activation sub-cortical and cortical investigated by MEG MEG or electroencephalography electroencephalography (EEG) can averaged event related field. In an EEG study (Porcaro et al. be found in Curio (2004). To assess fluctuations of SE-HFOs on a single trial (single evoked response) level is of notSE-HFOs possible averaged event related activation field. In anofEEG study (Porcaro et al. (2009)) the sequential sub-cortical and cortical be found in Curio (2004). To assess fluctuations SE-HFOs is demonstrated using sub-cortical the N20 response as be found in Curio (2004). To assess fluctuationsis of SE-HFOs the activation and on a to single (single evoked response) not possible due theirtrial weakness compared to otherlevel signals present in (2009)) (2009)) the issequential sequential activation of of sub-cortical and cortical cortical SE-HFOs the N20 response as on aa to single trial (single evoked level is possible constraint for andemonstrated ICA solution.using on single trial (single compared evoked response) response) level is not not possible SE-HFOs is demonstrated using the N20 response as due their weakness to other signals present in the MEG data. In thiscompared work a to statistical method for the SE-HFOs is andemonstrated constraint for ICA solution.using the N20 response as due to their weakness other signals present in due to theirdata. weakness compared to other signals present in constraint for an ICA solution. the MEG In this work a statistical method for the extraction of single trial SE-HFOs is successfully applied and constraint for anconcatenation ICA solution. temporal decorrelation (ecTD) the data. In this work method for the the epoch the MEG MEG of data. In trial this SE-HFOs work aa statistical statistical method for and the Here extraction single is successfully applied Here the epoch concatenation temporal decorrelation properties of the single trials are derived. extraction single trial SE-HFOs is successfully applied and ICA method (Zavala-Fernandez et al. (2012)) is used(ecTD) on the of the the epoch concatenation temporal decorrelation (ecTD) extraction single trialtrials SE-HFOs is successfully applied and Here properties of single are derived. Here the epoch concatenation temporal decorrelation (ecTD) ICA method (Zavala-Fernandez et al. (2012)) is used on the properties of the single trials are derived. SE-HFO signal section of the raw data. It consists of (Zavala-Fernandez et al. is on properties of the single trials are derived. The isolation of SE-HFO single trials is complicated by their ICA ICA method method (Zavala-Fernandez et raw al. (2012)) (2012)) is used used on the the SE-HFO signal section of the data. It consists of bandpass filtering the raw data between 450 and 1000 Hz, signal of the raw It of The single is complicated by their SE-HFO shortisolation durationof inSE-HFO comparison to trials a typical electrical stimulus SE-HFO signal section section ofdata thebetween raw data. data. Itandconsists consists of bandpass filtering the raw 450 1000 Hz, The isolation of SE-HFO single trials is complicated by their followed filtering by the concatenation of equal450 length windows The isolation ofinSE-HFO single trials is complicated by their bandpass the raw data between and 1000 Hz, short duration comparison to a typical electrical stimulus repetition rate in of comparison 3-9 Hz. ThetoSE-HFOs are easily stimulus seen in bandpass filtering the raw data between 450 and 1000 Hz, followed by the concatenation of equal length windows short duration a typical electrical covering the SE-HFO section of the data. The resulting time short duration in comparison toSE-HFOs a typical electrical stimulus by the of equal length windows repetition are easilyfor seen in followed only 5 msrate of of the3-9 111Hz. msThe interstimulus interval 9 Hz followed by SE-HFO the concatenation concatenation of data. equalThe length windows covering the of the resulting time repetition rate of 3-9 Hz. The SE-HFOs easily in series is subjected to section time-delayed (temporal) decorrelation repetition rate of 3-9 Hz. The SE-HFOs are are easilyforseen seen in covering the SE-HFO section of the data. The resulting time only 5 ms of the 111 ms interstimulus interval 9 Hz stimulation. This means that SE-HFOs are for sparsely covering the SE-HFO section of the (temporal) data. The resulting time series is subjected to time-delayed decorrelation only 5 ms of the 111 ms interstimulus interval 9 Hz (TDD)is(Belouchrani ettime-delayed al. (1997); (temporal) Ziehe et al.decorrelation (1998); de only 5 ms ofThis the 111 ms interstimulus interval for 9 Hz series subjected to stimulation. means that SE-HFOs are sparsely distributed in measured raw data corresponding to nonseries is subjected to time-delayed (temporal) decorrelation (TDD) (Belouchrani et al.TDD (1997); et al. (1998); stimulation. This SE-HFOs are Cheveigne et al. (2014)). was Ziehe successfully applied de in stimulation. This means meansrawthat that are sparsely sparsely (Belouchrani et (1997); Ziehe et de distributed dataSE-HFOs corresponding to nonstationarity: in To measured describe statistically the SE-HFO signal on a (TDD) (TDD) (Belouchrani et al. al.TDD (1997); Ziehe et al. al. (1998); (1998); Cheveigne et al. (2014)). was successfully applied de in distributed in measured raw data corresponding to nonthe study ofetN20m responses (Tang etsuccessfully al. (2002); Kishida et distributed in measured raw data the corresponding to nonstationarity: To describe statistically SE-HFO signal on a Cheveigne al. (2014)). TDD was applied 10 ms timeToscale a time dependent probability density Cheveigne etN20m al. (2014)). TDD was etsuccessfully applied in in the study of responses (Tang al. (2002); Kishida et stationarity: describe statistically the SE-HFO signal on a (2003); Tang et al. (2005); Sander et al. al. (2005)). stationarity: describe statistically the SE-HFO signal on a al. the study of N20m responses (Tang et (2002); Kishida 10 ms timeToscale a time dependent probability density the study ofTang N20m responses (Tang et (2002); Kishida et et al. (2003); et al. (2005); Sander et al. (2005)). 10 ms time scale a time dependent probability density 10 ms time scale a time dependent probability density al. (2003); Tang et al. (2005); Sander et al. (2005)). al. (2003); Tang et al. (2005); Sander et al. (2005)).
Copyright 2015 453 Hosting by Elsevier Ltd. All rights reserved. 2405-8963© 2015,IFAC IFAC (International Federation of Automatic Control) Copyright ©©2015 IFAC 453 Copyright © 2015 453Control. Peer review underIFAC responsibility of International Federation of Automatic Copyright © 2015 IFAC 453 10.1016/j.ifacol.2015.10.182
9th IFAC BMS 454 Tilmann H. Sander et al. / IFAC-PapersOnLine 48-20 (2015) 453–458 Aug. 31 - Sept. 2, 2015. Berlin, Germany
2. METHODS 2.1 Concatenated data and temporal decorrelation Given a superposition of n sources of the form ai si(t), where ai is a time independent field map of the m channel MEG system and si(t) describes the time dependence of the source, the measured data vector x(t) can be decomposed by a demixing matrix W: W x(t) = u(t).
(1)
Here the ui(t) are identical to the si(t) if disregarding an indeterminacy in scale and ordering. To estimate a W blindly ICA is frequently applied and n=m is assumed. One ICA method suitable for MEG data is the time-delayed decorrelation (TDD, TDSEP, SOBI) algorithm (Belouchrani et al. (1997); Ziehe et al. (1998); Congedo et al. (2008)), which is often classified as a blind source separation type algorithm as it relies solely on second-order statistics. In TDD a set of time-delayed covariance matrices is calculated, these matrices are symmetrised, and subsequently approximately diagonalised to obtain an estimate of W. As the SE-HFOs have a limited duration and a specific frequency the ecTD algorithm (Zavala-Fernandez et al. (2012)) derived from TDD seems suitable to extract them. Firstly the data are bandpass filtered between 450 and 1000 Hz. Then equal length windows covering the duration of the SE-HFO wavelet are concatenated to form a new continuous time series. Finally time-delayed covariance matrices of the concatenated time series are calculated and processed as in standard TDD.
Fig. 1. Illustration of input data for the temporal decorrelation ICA. a) TDD: Two time series x and y from a multi-channel measurement are used in the time-delayed covariance calculation C(x,y,τmax). b) Epoch concatenation TD: Concatenation of equal length time windows wi of x and y yields new time series xc and yc used for C(xc,yc,τmax).
To contrast ecTD with TDD the input data for the maximum delay τmax covariance matrix are illustrated in Fig. 1 for both algorithms. In Fig. 1a) two time series x and y of a multichannel measurement are sketched. The data range 0...T of x is correlated with the data range τmax...T+τmax of y (ranges indicated by red bar) to obtain the matrix element C(x,y,τmax) of TDD.
noise level. This avoids steps in the concatenated time series. For ecTD the τmax covariance matrix C(xc,yc,τmax) is calculated from concatenated time series such as xc and yc for the data ranges indicated in red. The SE-HFOs shown in Fig. 1b) are exaggerated and cannot be observed in real data.
Due to the wavelet like nature of SE-HFOs they can be modeled as s(t) = A(t) sin(2 π fHF t + φ)
(2)
The C(xc,yc,τk) matrix of the ecTD algorithm for M windows of length wlen is given by
with envelope A(t), frequency fHF, and phase φ. In this model f is about 600 Hz and the envelope A(t) has a support of less than 10 ms centered on the N20m peak. This model is motivated by physiological knowledge derived from electrical recordings, which show that the SE-HFOs are always associated with the N20m peak in a systematic fashion (Klostermann et al. (1999)).
Cij(xc,yc,τk) = Σ{l} xi(tl) xj(tl+ τk),
(3)
where the set of indices is {l} = w1,start,...,w1,start + wlen, w2,start,...,w2,start + wlen,wM,start,...,wM,start + wlen and contains M* wlen elements. The ecTD algorithm requires choosing three parameters: The window start wi,start relative to the stimulus, the window length wlen, and the maximum delay τmax. The window start and the window length are fairly obvious in the case of SEHFOs as the window should contain the complete SE-HFO wavelet. The length of the SE-HFO can be estimated from a conventional average. The choice of τmax is largely empirical,
The raw SE-HFO input data for ecTD are sketched in Fig. 1b) as a new concatenated time series xc is created from x using the data in the windows wi of length wlen indicated in blue. Any discontinuity between the windows will be obscured in the noise as the length of the windows extends beyond the SE-HFOs, i.e., the range where A(t) is above 454
9th IFAC BMS Aug. 31 - Sept. 2, 2015. Berlin, Germany Tilmann H. Sander et al. / IFAC-PapersOnLine 48-20 (2015) 453–458
455
an accepted rule of thumb is to set a desirable frequency resolution f and set τmax= 1/f.
but it was found that ecTD works best for 1/τmax >> wlen. Using a time-delay larger than the window length implies similar data and noise content in each window, i.e., stationarity. A rule of thumb for TDD states that the frequency resolution is 1/τmax (Sander et al. (2005)).
Once single trial SE-HFOs have been extracted they can be characterised using (2). Set values for fHF and a time range tstart ... tstop are chosen and then time independent average amplitudes Ai and phases φi are calculated for each trial (Delorme et al. (2003)). Only the three central oscillations of the wavelet are used to estimate Ai. Subsequently the single trials can be sorted by amplitude or phase or by a sequence of both. For visualization after sorting a single trial (ST) plot (Jung et al. (2001a); Tang et al. (2002)) will be used. In the ST plot multiple sections of a time series are displayed as horizontal scan lines stacked vertically with the function value coded as a grey shade between black and white.
Fig. 2. Conventionally averaged MEG result for an electrical stimulation at the left medianus nerve. The map sequence shows at 21 ms a dipolar pattern known as the N20m. Channel CH44 has the highest peak amplitude and small variations in the shape of the peak are the SE-HFOs directly visible in Fig. 3 after bandpass filtering. 2.2 Measurements and data processing In one healthy volunteer MEG was recorded using a 125channel whole head MEG system made by Eagle Technology (Yokogawa MEG Vision (1998)). To elicit a SEF N20m an electrical pulse of 100 μs duration and 9 Hz repetition rate was applied to the left medianus nerve with a strength of 12 mA selected to be above the motor threshold. The measurement lasted 670 s and the data were sampled at 5 kHz. All measurements were performed in accordance with the declaration of Helsinki. Data from other subjects showed similar results.
Fig. 3. Average over the highpass filtered data shown as a sequence of maps together with the result for a single channel (top part) and the HF-ecTD component map with the associated averaged time series (bottom part). The oscillations in both time series show the same phase and envelope and the SE-HFOs have a frequency of f ≈ 625 Hz.
Before the ecTD calculation the multi-channel raw data set was digitally filtered using a 450 to 1000 Hz bandpass filter with -10 dB damping points at 400 and 1100 Hz and no phase shifts in the passband. After digital filtering the epochs were concatenated to create a new data set (see section 2.1) using wstart = 16 ms post stimulus and wlen = 13 ms. This choice of wstart and wlen is valid for any other dataset of this type as the N20m occurs subject independent at 20 ms with a variation of at most 2 ms and therefore this choice covers the SE-HFOs riding on top of the N20m.
3. RESULTS 3.1 Extraction of SE-HFO signal The conventional N20m result is given in Fig. 2, which shows a map sequence from the average over 6000 epochs of recorded data and the time series of channel CH44, which has the highest peak amplitude at 20 ms. A 5 ms baseline starting at 5 ms after the stimulus was subtracted.
The set of time delayed covariance matrices was chosen empirically to be {τi} = {2, 4, 6,..., 2000 ms}, which requires to diagonalise 1000 matrices. The largest τi corresponds to a frequency of f = 1/τmax = 1/2000 ms = 0.5 Hz. A theoretical foundation for the choice of τmax for TDD is not known and
The map at t = 0 s is the stimulus artifact due to the electrical current pulse, at 16 ms no activity is seen, and a clear dipolar pattern is observed in the map at 21 ms. Later activity at 24 ms is more complex. The dipolar pattern at 21 ms is on the 455
9th IFAC BMS 456 Tilmann H. Sander et al. / IFAC-PapersOnLine 48-20 (2015) 453–458 Aug. 31 - Sept. 2, 2015. Berlin, Germany
right side of the map contralateral to the left stimulated side (orientation of map: top of map points to nose, left and right of map correspond to left and right hemisphere). A dipolar pattern as shown in the map at 21 ms was successfully approximated by an ECD source in the cortex for the same whole head MEG system (Sander et al. (2007)). To extract the SE-HFOs the raw data have to be filtered using the 450-1000 Hz bandpass mentioned previously (see section 2.2) and then averaged. The SE-HFOs appear then as a short wavelet around 21 ms in the average. This is shown in the top part of Fig. 3, where a map sequence and the single channel of the averaged filtered data are given. The SE-HFOs can be seen in Fig. 2 as small peak shape variations in the same channel (CH44) before filtering. The maps in Fig. 3 appear much noisier compared to the maps in Fig. 2, which is not surprising noting the different contour step for the maps in the Figs. 2 and 3 and the well known small amplitude of the SE-HFOs. The maps in Fig. 3 show a dipolar structure with a clear resemblance to the N20m map in Fig. 2. The time points of the maps were chosen to coincide with the oscillation extrema in the time series shown below the maps. This time series shows the channel Ch 44, which had the strongest oscillation amplitude. The location of this channel is indicated in the map. As the maps are aligned with the oscillation extrema the sign of the map reverses with a periodicity of 1.6 ms (= 625 Hz), the periodicity indicated for the time series. The first peak of the SE-HFOs appears at 17 ms, which is consistent with the literature (Curio et al. (1994); Hashimoto et al. (1996)).
Fig. 4. Comparison of power spectra of unaveraged data. Top: Bandpass filtered data of Ch 44. Middle: Concatenated bandpass filtered data of of Ch 44 (wlen = 13 ms). Bottom: The time series of the HF-ecTD component from Fig. 3. The peaks marked by red circles in the middle panel dominate the spectrum of the HF-ecTD component.
the concatenated data. The occurrence of several peaks in the spectrum of HF-ecTD can be explained as an amplitude modulation process. Assuming a single frequency f modulated by an envelope fm peaks at f+fm and f-fm result. The concatenated SE-HFOs can be viewed as such a type of signal with f ≈ 625 Hz and fM = 77 Hz. As the SE-HFOs do not have a simple sinusoidal envelope more peaks result in the spectrum. It was verified that changing the concatenation window length changes the peak separation. The ecTD algorithm apparently succeeded in isolating the SE-HFO signal due to its periodic signature above noise level in the concatenated signal. In agreement with that no component related to the SE-HFOs was obtained by applying TDD to the full length band pass filtered data using similar time-delay parameters as used in the ecTD calculation.
As an alternative method of extracting the SE-HFOs the ecTD algorithm was applied to the bandpass filtered data using the parameters described in the previous section. Only one of the 125 resulting component maps had a dipolar pattern, which is a sign of a brain source (Delorme et al. (2012)). This component map is shown in the bottom part of Fig. 3 and labeled as HF-ecTD. The average over its time series is similar to the Ch 44 time series as can be seen in the middle of Fig. 3 and therefore the component HF-ecTD is attributed to the SE-HFOs. All other ecTD components (not shown) have an irregular map typical of noise components. The HF-ecTD average time series covers only the range from 15 to 28 ms as this was the concatenation window for the ecTD algorithm. The HF-ecTD component map is similar to the map at 23 ms from the simple average (Fig. 3, top). The similarity of the two maps can be quantified by an angle α defined as cos α = a*b / ( |a|*|b| ) and α = 26o results, which means a high degree of similarity given the additional noise visible in the maps.
3.2 Single trial analysis The SE-HFOs were analysed subsequently using the model in (2). The parameters Ai and phase φi of each trial were estimated for fHF = 625 Hz and tstart ... tstop covering a 4.8 ms window around t = 23.6 ms post-stimulus. The choice of tstart and tstop was motivated by the envelope of CH44 in Fig. 2, which shows three peaks at the same amplitude followed by two peaks with a higher amplitude. The existence of two regimes in the SE-HFOs is well known (Klostermann et al. (1999); Haueisen et al. (2007)). This choice of tstart and tstop might be not suitable for other datasets and similarly the averaged SE-HFOs would need to be inspected.
An important property of the HF-ecTD component is presented in Fig. 4, which compares the power spectra of the unaveraged data for different processing stages. The spectrum of the bandpass filtered Ch 44 data does not have specific features, whereas four peaks separated by f = 1.0/ wlen = 77 Hz (= 1.0/13 ms) are observed above the background in the spectrum of the concatenated data (peaks indicated by red circles). Very distinct peaks are visible in the spectrum of HF-ecTD at the same positions as the peaks in
Then the trials were sorted by Ai (equation (2)), split into four groups, and sorted again within each group by phase φi. The result is shown in Fig. 5, where the ST plots (see section 2.2) of the high and low amplitude groups (top, middle) are shown together with the subaverages of the four groups 456
9th IFAC BMS Aug. 31 - Sept. 2, 2015. Berlin, Germany Tilmann H. Sander et al. / IFAC-PapersOnLine 48-20 (2015) 453–458
(bottom; orange = low amplitude, dark blue = high amplitude) for the time series of Ch 44 and of component HF-ecTD. For display purposes non-overlapping groups of ten trials were averaged and an 11 ms window is shown. The sign is changed deliberately in Fig. 5 between Ch 44 and HFecTD to facilitate a comparison of the ST plots.
457
Nevertheless the total average in Fig. 3 is very similar for the two ensembles of trials, i.e., both ensembles contain the same amount of signal power. The observation of a preferred phase in the HF-ecTD-ST plot is only possible through the separation of noise and signal in the ecTD calculation. Note that the four subgroup averages in Fig. 5 (bottom) of Ch 44 appear to have a higher noise level compared to those of HFecTD in agreement with a dominant noise removal effect. Only a single dipolar pattern, component HF-ecTD, was observed in the ecTD result despite the possibility to detect non-orthogonal sources using an ICA-type method. The oscillatory signal in the four SE-HFOs subaverages in Fig. 5 (bottom) is almost identical outside the sorting window, i.e., for t < 21 ms. Apparently the double sorting is relevant only for the 5 ms sorting window. This indicates that the weaker oscillations between 16 and 21 ms could originate from a different process compared to the stronger oscillations between 21 and 26 ms. This might be related to the early and late (before/after N20m) SE-HFO regimes observable in EEG and weaker in MEG data (Klostermann et al. (1999); Haueisen et al. (2007)). 4. CONCLUSION The successful extracting of an ensemble of SE-HFOs single trials by ecTD shows that ecTD is useful to extract transient oscillatory signals in multichannel recordings such as MEG. The ensemble extracted by ecTD is much better suited for further processing compared to the original ensemble as can be seen from the single trial analysis. Although this methodology was demonstrated on a single dataset it will be applicable to study group effects, because the SE-HFOs associated with the N20m can be observed in most subjects. The application of standard TDD (SOBI, TDSEP) to N20m data (Tang et al. (2002); Kishida et al. (2003); Tang et al. (2005); Sander et al. (2005)) is directly related to the present situation. The signal repetition and partitioning created through concatenation in ecTD is implicit in raw N20m data as the stimulation is periodic and the N20m neuronal response is strong enough enough to yield a distinct spectral pattern. For the N20m the temporal support has a sufficient length for temporal decorrelation to work, but for SE-HFOs the temporal extension is too small and concatenation is needed.
Fig. 5. Sorted ST plots for channel Ch 44 (left) and the time series of the HF-ecTD component (right). The trials were sorted by amplitude Ai and phase φi in four groups. Of the four groups the ST plots of the highest and lowest amplitude groups are shown (top, middle) together with the individual subaverages for the four groups (bottom).
Comparing upper and lower ST plot in Fig. 5 it is evident that the single trial SE-HFOs cover a wide range of amplitudes. Comparing the high amplitude ST plots of Ch 44 and HFecTD shows that a preferential phase exists in the single trials of HF-ecTD as most of the single trials have a peak aligned with the black vertical line added as a guide to the eye. The peaks of the HF-ecTD single trials form a sigmoid line, it is almost straight for Ch 44 corresponding to a uniformly distributed phase. The average over an oscillation with a uniformly distributed phase is zero for an infinite number of trials and mainly the amplitude distribution of the trials in Ch 44 leads to a non-zero average.
Clearly the known modulations of the SE-HFOs (Curio (2004)) should be investigated by extracting single trials using ecTD on a group of subjects, e.g., is the decay of the SE-HFO amplitude during sleep due to phase randomization or due to a decrease of amplitude? Newer studies (Fedele et al. (2015); Leiken et al. (2014)) have measured noninvasively electro-physiological brain signals around 1 kHz and ecTD might be suitable to get more insights into these data as well. ACKNOWLEGDMENTS This work is supported by BMBF grant 01 GO 0208 BNIC (Berlin Neuroimaging Centre) and DAAD scholarship A0421558.
Related to this are the grey scale bars on top of the Ch 44and HF-ecTD-ST plots showing that the single trials in Ch 44 have typically twice the amplitude compared to HF-ecTD. 457
9th IFAC BMS 458 Tilmann H. Sander et al. / IFAC-PapersOnLine 48-20 (2015) 453–458 Aug. 31 - Sept. 2, 2015. Berlin, Germany
REFRENCES
Kishida K., Kato K., Shinosaki K., Ukai S. (2003). Blind identification of SEF dynamics from MEG data by using decorrelation method of ICA. Proc. of the 4th Intl. Symposium on ICA nd BSS, Nara, 185-190. Klostermann F., Nolte G., Curio G. (1999). Multiple generators of 600 Hz wavelets in human SEP unmasked by varying stimulus rates. Neuroreport, 10, 1625-1629. Lee T.W. (1998). Independent Component Analysis - Theory and Applications. Kluwer Academic Publishers, Boston. Leiken K., Xiang J., Zhang F., Shi J., Tang L., Liu H., Wang X. (2014). Magnetoencephalography detection of highfrequency oscillations in the developing brain. Front. Hum. Neurosci., 8, 969. Makeig S., Debener S., Onton J., Delorme A. (2002). Mining event-related brain dynamics. TRENDS in Cogn. Sciences, 8, 204-210. Porcaro C., Coppola G., DiLorenzo G., Zappasodi F., Siracusano A., Pierelli F., Rossini P.M., Tecchio F., Seri S. (2009). Hand somatosensory subcortical and cortical sources assessed by functional source separation: An EEG study. Hum. Brain Mapp., 30, 660-674. Roberts S., Everson R. (eds.) (2001). Independent Component Analysis - Principles and Practice. Cambridge University Press, Cambridge. Sander T.H., Burghoff M., Curio G., Trahms L. (2005). Single evoked somatosensory MEG responses extracted by time delayed decorrelation. IEEE Trans. on Signal Processing, 53, 3384-3392. Sander T.H, Burghoff M., Van Leeuwen P., Trahms L. (2007) Application of decorrelation independent component analysis to biomagnetic multi-channel measurements. Biomedizinische Technik 52: 130–136. Stone J.V. (2002). Independent component analysis: An introduction. TRENDS in Cogn. Sciences, 6, 59-64. Tang A.C., Pearlmutter B.A., Malaszenko N.A., Phung D.B. (2002). Independent components of magnetoencephalography: Single-trial response onset times. Neuroimage, 17, 1773-1789. Tang A.C., Sutherland M.T., McKinney C.J. (2005). Validation of SOBI components from high-density EEG. Neuroimage, 25, 539-553. Yokogawa MEG Vision (1998). Distributed by YOKOGAWA / Japan. http://www.yokogawa.com/us/technical-library/whitepapers/megvision-magnetoencephalograph-system-andits-applications.htm (accessed 23.1.2015). Zavala-Fernandez H., Orglmeister R., Trahms L., Sander T.H. (2012). Identification enhancement of auditory evoked potentials in EEG by epoch concatenation and temporal decorrelation. Comput Methods Programs Biomed. 108, 1097-105. Ziehe A., Müller K.-R. (1998). TDSEP - an efficient algorithm for blind separation using time structure. In L. Niklasson et al. (eds.), Proceedings of the 8th ICANN, 675-680. Springer.
Belouchrani A., Abed-Meraim K., Cardoso J.-F., Moulines E. (1997). A blind source separation technique based on second-order statistics. IEEE Trans. Sig. Proc., 45, 434444. Choi S., Cichocki A., Park H.-M., Lee S.-Y. (2005). Blind source separation and independent component analysis: A review. Neural Inf. Proc., 6, 1-57. Cichocki A., S.Amari S. (2002). Adaptive Blind Signal and Image Processing . John Wiley and Sons, New York. Congedo M., Gouy-Pailler C., Jutten C. (2008). On the blind source separation of human electroencephalogram by approximate joint diagonalization of second order statistics. Clin. Neurophysiol., 119, 2677-2686. Curio G., Mackert B.-M., Burghoff M., Koetitz R., AbrahamFuchs K., Härer W. (1994). Localization of evoked neuromagnetic 600Hz activity in the cerebral somatosensory system, Elec. Clin. Neurophysiol., 91, 483-487. Curio G. (2004). Ultrafast EEG Activities. In E. Niedermeyer and L. da Silva (eds.), Electroencephalography - Basic Principles, Clinical Applications, and Related Fields, 495-504. Lippincott, Williams and Wilkins, Phialdelphia. de Cheveigné A., Parra L.C. (2014). Joint decorrelation, a versatile tool for multichannel data analysis. Neuroimage, 98, 487-505. Delorme A., Makeig S. (2003). EEGLAB: An open source toolbox for analysis of single-trial EEG dynamics including independent component analysis. J. of Neurosc. Meth., 134, 9-21. Delorme A., Palmer J., Onton J., Oostenveld R., Makeig S. (2012). Independent EEG sources are dipolar. PLoS One, 7, e30135. Fedele T, Scheer HJ, Burghoff M, Curio G, Körber R. (2015). Ultra-low-noise EEG/MEG systems enable bimodal non-invasive detection of spike-like human somatosensory evoked responses at 1 kHz. Physiol. Meas., 36, 357-368. Hashimoto I., Mashiko T., Imada T. (1996). Somatic evoked high-frequency magnetic oscillations reflect activity of inhibitory interneurons in the human somatosensory cortex. Elec. Clin. Neurophysiol., 100, 189-203. Haueisen J., Leistritz L., Süsse T., Curio G., H.Witte H. (2007). Identifying mutual information transfer in the brain with differential-algebraic modeling: Evidence for fast oscillatory coupling between cortical somatosensory areas 3b and 1. Neuroimage, 37, 130-136. Haykin S. (ed.) (2000). Unsupervised Adaptive Filtering. Volume I: Blind Source Separation. John Wiley and Sons, New York. Hyvärinen A., Karhunen J., Oja E. (2001). Independent Component Analysis. John Wiley and Sons, New York. James C.J., Hesse C.W. (2005). Independent component analysis for biomedical signals. Physiol. Meas., 26, R15R39. Jung T.P., Makeig S., McKeown M.J., Bell A.J., Lee T.W., Sejnowski T. (2001). Imaging brain dynamics using independent component analysis. Proc. of IEEE, 89, 1107-1122.
458