Single-Trial Estimation and Analysis of PVEP Based on Independent ComPonent Analysis

TSINGHUASCIENCE AND TECHNOLOGY;-508ISSN 1007-0214 21/21pp503。 508Volume 6. Number 5, December 2001Single-Trial Estimation and Analysis of PVEP Based onIndependent Component Analysis'HONG Bo (洪波)， YANG Fusheng (杨福生), PAN Yingfu (潘映辐)'TN ATANG Qingyu (唐庆玉),CHENKui(陈奎)| ，TIE Yanmei (铁艳梅)'Department of Electrical Enginering，Tsinghua University, Beijing 10084. China;十Beijing Friendship Hospital, Beijing 100050, ChinaAbstract: A novel method based on the extended intomax approach of independent component analysis(ICA) is proposed for the single trial estimation of multi-channel pattern-reversal visual evoked potential(PVEP>. The ctinical aplications Hlustrate the encouraging performance of the method. The number of trialsneeded is reduced to 3, but the estimated VEP complex is clearer than that obtained by conventional ensembleaveraging with 50 trials. By analyzing the time course and spatial pattern of the independent components (ICs)of the multi channel PVEP, one component is found to be reasonably related to the P100 peak.Key words : independent component analysis; patternr reversal visual evoked potential (PVEP); single trialestimation; P100 peak1techniques including adaptive filteringl), waveletIntroductiondenoising-i, and neural network filteringPattern. reversal visual evoked potential (PVEP) is .Recently，blind source separation bywidely used in clinical diagnosis and surgery ，independent component analysis ( ICA ) hasbecause it can rflect the result of the visualreceived much attention because of its potentialprocess consisting of light perception at theapplications in biomedical signal processing such asretinal, transmittance along the visual pathways,electrocardiograph ( ECG )[5]，electroencephalo-and the final cortical presentation. In practice, thegraph (EEG )C5-8] and functional MRI ({MRI)[9].PVEP is a transient signal contaminated by theThe goal of ICA is to recover independent sourcesongoing background electroencephalograph (EEG)from observations that are unknown linearand other noises with poor signal to noise ratio. Tomixtures of the unobserved independent sourceobtain a clear PVEP complex, the ensemblesignalsC1o]. Neurophysiological studies have shownaveraging method is commonly used with a largethat the evoked potential is generated innumber of trials. However, the PVEP signal canconjunction with synchronous activity in pyramidalnot be considered as a deterministic process sincecells in the activated cortical area (the source) andits components are known to vary both inthe volume conduction through the cerebrospinalamplitude and latency from stimulation tofluid，skull and scalp is thought to be linear.stimulation'', so the ensemble averaging approachTherefore, ICA can be utilized to get a better viewmay lead to the loss of dynamic information. Inof the electric activities in the brain by demixingaddition, the number of trials should be reduced inthe measured scalp EEG into its statisticallyclinical situations to avoid fatigueandindependent components. More precisely, ICAaccommodation of the patient's visual system. Fordecomposes the PVEP recorded by multiplethese reasons , many efforts have been contributedelectrodes on the scalp into a sum of componentsto the development of single-trial estimationwith fixed spatial distributions and sparselyactivated, maximally independent time courses ，,Received: 2000 0620without attempting to directly specify where in the* Supporied by the National Natural Science Foundationbrain these activities were generated.of China (No. 59937160)based method中国煤化工"fYRCNMHGTsinghua Science and Technology, December 2001, 6(5); 503 - 508504was first applied to the enhancerent and analysisstatistical higher order moments and try toof multi -channel PVEPs. A novel method based oneliminate them, which is called the cumulant-basedthe optimal selection of independent componentsmethod. The other is to apply the nonliner(ICs) was utilized to get a clear estimate of the squashing function to produce the necessasry higherPVEP complex in a single trial. In the clinicalorder statistics through its Taylor seriesapplication, the proposed method was shown toexpansion.Consideringthecomputationalgive better performance than the ensemble-efficiency, the later approach is more convenient toaveraging method. In addition, the single-trialapply than the former one. In most cases, theapproach of PVEP estimation allowed us toneural network learning approach is used toinvestigate the dynamic change of NPN complex ofoptimize the objective function that implies theVEP from trial to trial. The spatial patternindependence between signals. In this frameworkembedded in the demixing matrix derived by ICA(Fig. 2), each output component u,is fed to amay also help us determine the ICs that havenonlinear function g; ( ) which was used topossible physiological meaning.approximate the calculation of the higher ordercumulant. It can be proved that the ideal form of1 Methodg;( ) is the cumulative density function (c.d.f.) ofthe distributions of the independent sources syo1.1 Independent component analysis.However, many authors have found in their workThe ICA idea is closely related with the problem ofthat the choice of the nonlinear function gj( ) isblind source separation (BSS). A vector ofnot very critical1.121. For example, the Sigmoid orindependentsourcess(t)=[(,"..sw(t)]TTanh functions may serve as good candidates forpropagates through a medium and mixtures ofg,( ) when the sources have super-Gaussianthem are picked up by N sensors. The observeddistribution.vector is denoted as x=As，where A is a squarematrix. The BSS problem is to identify the matrix，+AA, or to find a matrix W such that u=Wx=WAs isan estimate of the source signal 8 (the entries of8( )vector u may be permuted and rescaled versions ofAdaptivethe entries of s) (Fig. 1). Since the sources s,(t)，leamingj= 1,2..N are independent of each other, thejoint probability density function (p.d.f.) ofs isFig. 2 Neural network learning approach to achieve theindependence of output signalsthe product of the marginal p. d. f. of eachcomponent :In the framework ilustrated in Fig. 2, Bellf,- if,(1)and Sejnowski used the sigmoid function to derivethe stochastic gradient learning rules fomaximization of the mutual information betweenthe input x and the nonlinear output y of a neural20 \A/ O\w/ Onetwork, which leads to the independence betweenthe linear outputs u13]. The learning rule can be30expressed as :oW =∈[(W")-'+ (I - 2y)x](2)nOwhere E is the learning rate. This informationFig. 1 Blind source separation ( BSS ) problem:theoretic approach is commonly called infomax ICAblindly separates the mixture (x) of the independentalgorithm. By applying the concept of naturalsource (s)gradientla2], the learning rule can be modified as:OW =e[(WT)-'+ (1- 2y)x"]w"W =Independence means the absence of all highere[1 + (1 - 2y)u"]W(3)orderstatisticalmomentsbetweenSince matrix inversion is avoided in Eq. (3), it hascomponents. The requirement is stricter than thatbetter computational performance and convergenceof the principle component analysis (PCA) inspeed.which only second-order correlation is removed.However, the algorithm can only be used toThere are two ways to achieve the independenceseparateourceswithsuper-Gaussianbetween signals. One is to evaluate all the中国煤化工y of extending theTHCNMHG505HONGeo (洪波) et al; Single Trial Esimation and Analysis of PVEP ".learning rule to sources with either sub- or super-waveforms) which are related with PVEP (foroftheexample, those ICs which show prominent peaksGaussian distributions is to use the sign.normalized kurtosis of the output uas the flag foror valleys during the time period from 0 to 300 ms)switching between these two conditionsho.and set the remaining ICs to zero. By so doing, .The algorithm is commonly referred as themost parts of the uncorrelated noises as well as thespontaneous EEG will be eliminated. If the newextended infomax ICA.vector so obtained is denoted as u', then the1.2 Single-trial estimation based on the extendedPVEP component in the scalp records can binfomax ICAreconstructed by multiplying the new vector uWhen the infomax ICA framework is applied towith W-':x'(t)= W-'u' (t)multi-channel PVEP data, the observed vector x isformed by the time series recorded at the differentAs expected, in the reconstructed data.EEG eletrodes. We retrieve the unknown mixinglear NPN complexes are found at occipitalmatrix A by estimating its inverse matrix W usingchannels, such as 01, O2 and Oz.the extended infomax ICA. The wholeThe ICA provides not only the time course butcomputation consists of two stages: the learningalso the spatial distribution of the ICs. Thestage to find the optimal matrix W and the workingcolumns of the inverse matrix W-' give thestage to enhance the PVEP by_ choosing theprojection strengths of the respectiveICAindependent cormponents related to PVEP.components onto the scalp electrodes. The reasonIn the learning stage, the observed data x isis explained as follows:divided into blocks and fed into the neural networkIf a single ICA component u, is chosen from u'consecutively and cyclically until matrix Wnd all the other components are set to zero,converges. Before that, a prewhittening step isEq. (5) can be rewritten as (where rui;jis thei, j-thused to improve the convergence speed and theentry of W-1):computational stability, in which mean values are「x(I)「W1:Wij : wIw] 「eliminated from each record of the training data:and then processed by a decorrelation technique tox;(t) =| WnVsW;N| u,(t)transform their covariance matrix into an identity|matrix. The data so obtained is then sent to theLxr(I)]WUMneura! network to find the optimal W matrix.(6)Meanwhile, sorme commonly used techniques fori.e., x(t) = Ww,u,(I)，i-l- Nneural network training, such as momentum termswhere z(t) is the projection of the ICA componentand adaptive learning rate adjustment , are properlyuj(t) onto the i-th electrode and wi,is its projectionintegrated into the learning scheme.strength. Thus, the column vector w, = [Wyj, W2,In the working stage (Fig. 3), we get the“,WwjJ is 1he projection strength over all the scalpoutput vector u by multiplying the data vector xelectrodes which can be used to reveal the spatialwith the matrix W found from the learning stage :pattern of u, (t) (at any time instant l) byu= Wx(4)topographic interpolation. The spatial pattern doesnot change with time. The temporal and spatialinformation of independent components of multi-二rmodificationw'二channe! EEG data can be combined to understandtheir physiological meaning[6^8].Fig. 3 Working stage diagram2 ResultsEach row of u is theactivation waveform of aCurrently,traditionalensembleaveragingcertain independent component (IC). But, it musttechniques are widely used in clinical applicationsbe emphasized that the PVEP we need ( whichto extract PVEP with 100200 trials usuallyexists in the scalp EEG records ) cannot be pickedneeded for the averaging. Even in some high-endup directly from the component of the vector u.instruments with denoising techniques, 50 trialsSome components of u are PVEP sources, but theyre the minimum requirement for an acceptableare located in the brain. To obtain the PVEP in thePVEP extraction，Besides the smearing of usefulscalp EEG records, further steps must be taken.information and the patient's uncomfortable feelingHere we choose those components of u (from theircaused中国煤化工extractionYHCNMHGTsinghua Science and Technology, December 2001, 6(5): 503 -508506of PVEP by ensemble averaging is unable to get( EBNeuro Co.). And the signal processingthe dynamic change of the PVEP waveform withdescribed here was performed on a PC. .the progression of stimulation which is often2. 2 Result on normal subjects and patientsdesired in neurosurgery and neurophysiologystudies. Our single-trial estimation and analysisOur single-trial estimation method based on ICAmethod based on ICA shows the potential ofwas tested using the acquired 16 -channel PVEPremoving these shortcomings.data. Figure 4 ilustrates the estimation result forone normal case: the left column is the scalp2. 1 Data acquisitlonpotentials o[ O1-Fpz, Oz Fpz and O2-Fpz (fromThe study was performed on a group of 72 patientstop to bottom)，the middle column is the ensemblefor clinical PVEP evaluation (57 healthy and 15averaging result (with 50 trials) and the rightabnormal according to catamnesis). The left andcolumnn is the result of our 1CA method.right hemifield checkerboard pattern reversal( Although our method is aimed at single-trialevoked potential was recorded by 16 Ag-AgClestimation，in the current study, we had toelectrodes placed according to the internationalperform all of our analysis with 3-trial averaged10/20 system. Linked ears served as the reference.data, which is the lowest limit of averaging in ourThe traditional black and white checker boardEEG acquisition system). As shown in Fig.4, thepattern reversed at a rate of 2 Hz.ICA based method gives clearer NPN complex thanPVEP was digitized at 1024 Hz, after which thethe ensemble averaging method, so that theremarkable eye movement and technical artifactsamplitude and the latency of each component of thewere removed. All the data acquisition wasNPN complex can be easily measured for clinicalperformed on a GALELIO EEG/VEP systemdiagnosis.Recorded PVEPEnsemble averagingICA enhancing言swwwwWwwM-20100 2003000 20030011 mst1 msFig. 4 ICA enhancing result comparing with the ensemble avereging methodThe P100 latency is an important indexwas shown to be abnormal by both methodsbut itcommonly used in clinical diagnosis. The P100was easier to quantify the abnormal latency fromlatencies of the 57 normal subjects were examinedthe ICA based result. Another case was a cerebralby both the traditional ensemble averaging methodinfarction patient whose P100 component wasand our ICA based method, The results show goodmuch delayed and with much lower amplitude thanconsistency between them: the ensemble averagingnormal. As shown in Fig. 6, the averaging resultresult is (106. 3士4. 5) ms and that of the ICA(left column) does not clearly show the P100based method is (106. 8土5.1) ms.which is clear in the ICA based result ( rightClinically, patients with P100 latency longercolumn).than the normal value plus three times theTemporal-spatial pattern of P100standard deviation are considered as abnormal. Forexample, the P100 latency of multiple cerebralmentioned above, both thtemporalsclerosis patients is usually found to be muchinformation and the spatial information of thelonger than the normal value. Figure 5 shows theindependent component can be obtained by ICA.PVEP of a multiple cerebral sclerosis patientOur current study investigates the physiologicalestimated by averaging (left column) and ICAmeaning of the independent component combiningbased (right column) methods. The P100 latency中国煤化工-rns. As shown in .*TYHCNMHGHONGBeo (洪波)et al: Sigle-Tial Estimation and Analysis of PVEP -..507Ensemble averagingICA enhancingO1-Fpz}O1-FpzOz-Fpz140ms/↓02-FpzW/02-Fpz'10000302003001/mst/. msFig.5 Delayed P100 in a multple sclerosis patient's PVEP is revealed clearly by the ICA enhancingOI-Fpz|主10-(几-1010158 ms! QOz-FPpzI 02-Fpz.Q2-FpzFig.6 Abnormal NPN complex in a cerebral infarection patient's PVEP revealed by ICA eahancingFig.7，the waveform shown at the top is anspatial pattern may be closely related with P100. Aindependent component that shows clear activalionrelation between the time course of the ICAaround 100 ms.Moreover, this specificcomponent and its spatial pattern is thusindependent component can always be found in allestablished.of our normal cases. If the corresponding column2.4 Dynamic variation of PVEP waveformof the matrix W 1 is mapped onto the scalp bymeans of EEG topography, a simple spatial patternThe ensemble averaging method discards thewith the highest level of activation on the occipitaluncorrelated noises as well as the variance amonglobe will appear which implies that this specificthe trials, so the dynamic variations of the VEPfrom trial to trial are lost. If the single-trialestimationbecomes possible, such dynamicinformation can then be captured. As a roughIC2indication of the ability of the 1CA method tocapture the dynamic variation, we ilustrate thechange here by comparing the single -trial PVEPestimation before and after 100 cycles o[ thecheckerboard pattern-reversal stimulation ( Fig.8). The after-stimulation PVEP (solid line) has alonger P100 latency and a wider NPN complexthan before the stimulation (dotted line ), whichIC2 Toporapindicates that the long-term repeating stimulationFig.7 Plausible P100 activation in an independentmay change the shape o[ the underlying PVEP.component and its spathal distribution on the scalpThis change. can be considered as a sign of the中国煤化工MYHCNMHGTsinghua sciencee and Techmology, December 2001, 6(5): 503 - 508508fatigue and the accommodation of the visual[2] Thakor N V. Adaptive filtering of evoked potentials.IEEE Trans on BME, 1987, BME -34 (1); 6-12.system.[3] Hong Bo, Yang Fusheng. Tang Qingyu. Single-trial10estimation of visual evoked potential combiningOz-Fpzwavelet shrinkage and adaptive neural network filter.n: Proceedings Satelite Symposium o[ the 20thAnnual International Conference of the IEEE EMBS,1998，129-132.S[4]FungKsM，LamFK,ChanFHY,etal.Adaptive neural network filter for visual evoked-stpotential estimation. In: Proceedings IEEEInternational Conference on Neural Networks, 1995,-105200005; 2293 -296.11 ms[5] Cardoso J F. Multidimensional independentFig.8 Variation of PVEP waveform before (ottedcomponent analysis. In: Proceedings IEEEline) and after (solid line) 100-cycle stimulationInternationa! Conference on Acoustics, Speech andSignal Processing, 1998, IV: 1941-1944.3 Conclusions and Discussion[6] Jung T P, Makeig S, Humphris C, ct al. Removingeletroencephalographic artifacts by blind sourceIn conclusion， a single-trial PVEP es timationseparation. Psychophysiology, 2000, 37; 163-178.method based on ICA was proposed and tested[7] Makeig S, Jung T P, Bell A J， et al. Blindhere. The simulation and clinical results indicateseparation of auditory event-related brain responsesthe potential application of ICA in the temporal-into independent components. Proc Natl Acad Sci,spatial analysisof the evoked potential. TheUSA, 1997, 94: 10 979- 10 984.overall framework of the application of ICA to[8] Makeig s，Westerfield M，Jung T P, et al.EEG analysis may be briefly summarized as: letFunctionally independent components of the latethe rows of the original data matrix x be the EEGpositive event related potential during visual spatialsignals recorded at different electrodes, the rowsattention. Journal of Neuroscience, 1999, 19 (7);of the ICA output matrix u= Wx will then give the2665 - 2680.time courses of activation of the decomposed[9] McKeown M J, Jung T P, Mekeig s, et al. Spatiallyindependent components while the columns ofindependent activity patterns in functional magneticmatrix W-1 give the projection strengths of theresonance imaging data during the Stroop color-respective independent components onto the scalpnaming task. Proc Natl Acad Sci, USA, 1998, 95:803 - 810.electrodes.ICA is a newly developed signal processing[10] Comon P. Independent component analysis- -A newconcept? Signal Processing, 1994, 36: 287- 314.technique that still has many open questions for[11] Cichocki A， Moszczynski L. Robust learninginvestigation. At present, in the ICA theory, thalgorithm for blind separation of signals. Electronicsmixing should be linear and instantaneous (noLetters, 1994, 30(17); 1386 - 1387propagation delays and no convolution) and the[12] Cardoso J F. Blind signal separation; statisticalnumber of sensors (i. e.，electrodes in EEG case)principles. Proc IEEE, 1998, 86(10): 2009 - 2025.should not be less than the oumber of sourcesle.183.[13] Bell A J， Sejnowski T J， An infornationThese limitations may not be completely satisfiedmaximization approach to blind separation and blindin actual EEG analysis. Although numericaldeconvolution. Neural Computation, 1995， 7 (6):simulations and experimental results have1129 - 1159.confirmed that we can capture the time course and[14] Lee T W. Girolami M, Sejnowski T J. Independentscalp topography of the temporally independentcomponent analysis using an extended infomaxsources especially in the case of evoked potentialalgorithm for mixed subgaussian and supergaussiananalysis[7.8], the limitations mentioned above arsources. Neural Computation， 1999， 11 (2 ):still worthy of further study.409- 433.References[1] McGillem C D, AunonJI, YuK B. Signal and noisein evoked brain potential. IEEE Trans on BME. .1985, BME-32(12): 1012- 1016.中国煤化工MYHCNMHG