Determining patterns in neural activity for reaching movements using nonnegative matrix factorization

MISSING IMAGE

Material Information

Title:
Determining patterns in neural activity for reaching movements using nonnegative matrix factorization
Series Title:
EURASIP Journal on Applied Signal Processing
Physical Description:
Book
Language:
English
Creator:
Kim, Sung-Phil
Rao, Yadunandana N.
Erdogmus, Deniz
Sanchez, Justin C.
Nicolelis, Miguel A. L.
Principe, Jose C.
Publisher:
BioMed Central
Hindawi Publishing Corporation
Publication Date:

Notes

Abstract:
We propose the use of nonnegative matrix factorization (NMF) as a model-independent methodology to analyze neural activity. We demonstrate that, using this technique, it is possible to identify local spatiotemporal patterns of neural activity in the form of sparse basis vectors. In addition, the sparseness of these bases can help infer correlations between cortical firing patterns and behavior.We demonstrate the utility of this approach using neural recordings collected in a brain-machine interface (BMI) setting. The results indicate that, using the NMF analysis, it is possible to improve the performance of BMI models through appropriate pruning of inputs.

Record Information

Source Institution:
University of Florida
Holding Location:
University of Florida
Rights Management:
All rights reserved by the source institution.
System ID:
AA00007613:00001


This item is only available as the following downloads:


Full Text

PAGE 1

EURASIPJournalonAppliedSignalProcessing2005:19,3113c 2005HindawiPublishingCorporationDeterminingPatternsinNeuralActivityforReachingMovementsUsingNonnegativeMatrixFactorizationSung-PhilKimDepartmentofElectricalandComputerEngineering,UniversityofFlorida,Gainesville,FL32611,USAEmail:phil@cnel.u.eduYadunandanaN.RaoMotorolaInc.,FL,USAEmail:yadu@cnel.u.eduDenizErdogmusDepartmentofComputerScienceandBiomedicalEngineering,OregonHealth&ScienceUniversity,Beaverton,OR97006,USAEmail:derdogmus@ieee.orgJustinC.SanchezDepartmentofPediatrics,DivisionofNeurology,UniversityofFlorida,Gainesville,FL32611,USAEmail:justin@cnel.u.eduMiguelA.L.NicolelisDepartmentofNeurobiology,CenterforNeuroengineering,DukeUniversity,Durham,NC27710,USAEmails:derdogmus@ieee.org;nicoleli@neuro.duke.eduJoseC.PrincipeDepartmentofElectricalandComputerEngineering,UniversityofFlorida,Gainesville,FL32611,USAEmail:principe@cnel.u.eduReceived31January2004;Revised23March2005WeproposetheuseofnonnegativematrixfactorizationNMF)asamodel-independentmethodologytoanalyzeneuralactivity.Wedemonstratethat,usingthistechnique,itispossibletoidentifylocalspatiotemporalpatternsofneuralactivityintheformofsparsebasisvectors.Inaddition,thesparsenessofthesebasescanhelpinfercorrelationsbetweencorticalringpatternsandbehavior.Wedemonstratetheutilityofthisapproachusingneuralrecordingscollectedinabrain-machineinterface(BMI)setting.Theresultsindicatethat,usingtheNMFanalysis,itispossibletoimprovetheperformanceofBMImodelsthroughappropriatepruningofinputs.Keywordsandphrases:brain-machineinterfaces,nonnegativematrixfactorization,spatiotemporalpatterns,neuralringactiv-ity.1.INTRODUCTIONBrain-machineinterfacesBMIsareanemergingeldthataimsatdirectlytransferringthesubjectsintentofmove-menttoanexternalmachine.Ourgoalistoengineerde-vicesthatareabletointerpretneuralactivityoriginatinginthemotorcortexandgenerateaccuratepredictionsofhandposition.IntheBMIexperimentalparadigm,hundredsofmicroelectrodesareimplantedinthepremotor,motor,andposteriorparietalareasandthecorrespondingneuralactiv-ityisrecordedsynchronouslywithbehavior(handreach-ingandgraspingmovements).Spikedetectionandsortingalgorithmsareusedtodeterminetheingtimesofsin-gleneurons.Typically,thespike-timeinformationissum-marizedintobincountsusingshortwindows100millisec-ondsinthispaper).Anumberoflaboratoriesincludingourownhavedemonstratedthatlinearandnonlinearadaptivesystemidenticationapproachesusingthebincountinput

PAGE 2

3114EURASIPJournalonAppliedSignalProcessing canleadtoBMIsthate ectivelypredictthehandpositionandgraspingforceofprimatesfordi erentmovementtasks[ 1 2 3 4 5 6 7 8 ].Theadaptivemethodsstudiedthusfarincludemovingaveragemodels,time-delayneuralnetworksTDNNs),Kalmanerandextensions,recursivemultilayerperceptrons(RMLPs),andmixtureoflinearexpertsgatedbyhiddenMarkovmodels(HMMs).BMIsopenupanimportantavenuetostudythespatio-temporalorganizationofspiketrainsandtheirrelationshipswithbehavior.Recently,ourlaboratoryhasinvestigatedthesensitivityofneuronsandcorticalareasbasedontheirroleinthemappinglearnedbytheRMLPandtheWienerter[7 ].Weexaminedhoweachneuroncontributestotheoutputofthemodels,andfoundconsistentrelationshipsbe-tweencorticalregionsandsegmentsofthehandtrajectoryinareachingmovement.Thisanalysisindicatedthat,dur-ingeachreachingaction,specicneuronsfromtheposteriorparietal,thepremotordorsal,andtheprimarymotorregionssequentiallybecamedominantincontrollingtheoutputofthemodels.However,thisapproachreliesondeterminingasuitablemodel,becauseitexplicitlyusesthelearnedmodeltoinferthedependencies.Inthispaper,weproposeamodel-independentmethod-ologytostudyspatiotemporalpatternsbetweenneuronalspikesandbehaviorutilizingnonnegativematrixfactoriza-tion(NMF)[9 10 ].Initsoriginalapplications,NMFwasmainlyusedtoprovideanalternativemethodfordetermin-ingsparserepresentationsofimagestoimproverecognitionperformance[10 11 ].AvellaandTreschhavealsopro-posedanextensionofNMFtoextracttime-varyingmusclesynergiesfortheanalysisofbehaviorpatternsofafrog[12 ]. ThenonnegativityconstraintsinNMFresultintheunsuper-visedselectionofsparsebasesthatcanbelinearlycombined(encoded)toreconstructtheoriginaldata.OurhypothesisisthatNMFcansimilarlyyieldsparsebasesforanalyzingneu-ralringactivity,becauseoftheintrinsicnonnegativityofthebincountsandthesparsenessofspiketrains.TheapplicationofNMFtoextractlocalfeaturesofneu-ralspikecountsfollowsthemethodofobtainingsparsebasestodescribethelocalfeaturesoffaceimages.Thebasisvec-torsprovidedbyNMFandtheirtemporalencodingpatternsareexaminedtodeterminehowtheactivitiesofspecicneu-ronslocalizetoeachsegmentofthereachingtrajectory.Wewillshowthattheresultsfromthismodel-independentanal-ysisoftheneuronalactivityareconsistentwiththepreviousobservationsfromthemodel-basedanalysis.2.NONNEGATIVEMATRIXFACTORIZATIONNMFisaproceduretodecomposeanonnegativedatamatrixintotheproductoftwononnegativematrices:basesanden-codingcoe cients.Thenonnegativityconstraintleadstoaparts-basedrepresentation,sinceonlyadditive,notsubtrac-tive,combinationsofthebasesareallowed.Ann m nonnegativedatamatrixX ,whereeachcolumnisasamplevec-tor,canbeapproximatedbyNMFasX = WH + E ,(1whereE istheerrorandW and H havedimensionsn r and r m ,respectively.W consistsofasetofr basisvectors,whileeachcolumnofH containstheencodingcoe cients foreverybasisforthecorrespondingsample.Thenumberofbasesisselectedtosatisfyr n + m
PAGE 3

DeterminationofNeuralFiringPatternsUsingNMF3115 Table1:Distributionofneuronsandcorticalregions. Monkey-1 Monkey-2 Regions PPM1(area1PMdM1(area2 M1PMdNeurons 1 3334 5455 8182 104 1 3738 54 3.1.DatapreparationSynchronous,multichannelneuronalspiketrainswerecol-lectedatDukeUniversityusingtwofemaleowlmonkeys(Aotustrivirgatus):Belle(monkey-1)andCarmenmonkey-2). 1 Microwireelectrodeswereimplantedincorticalregionswheremotorassociationsareknown[1 13 ].Duringtheneuralrecordingprocess,uptosixty-fourelectrodeswereimplantedinposteriorparietal(PP)-area1,primarymotor(M1)-area2,area4,andpremotordorsal(PMd)-area3,eachreceivingsixteenelectrodes.Fromeachelectrode,onetofourneuronscanbediscriminated.Theringtimesofindividualneuronsweredeterminedusingspikedetectionandsortingalgorithms[14 ]andwererecordedwhiletheprimateper-formeda3Dreachingtaskthatconsistsofareachtofoodfollowedbyeating.Theprimateshandpositionwasalsorecordedusingmultipleberopticsensorswithasharedtimeclock)anddigitizedwitha200Hzsamplingrate[1 ]. Thesesensorswerecontainedintheplasticstripofwhichbendingandtwistingmodiedthetransmissionofthelightthroughthesensorsinordertorecordpositionsin3Dspacemoreaccurately.Theneuronalingtimeswerebinnedinnonoverlappingwindowsof100milliseconds,representingthelocalringrateforeachneuron.Inthisrecordingsessionofapproximately20minutes(12000bins),104neuronsformonkey-1and54neuronsformonkey-2couldbediscrimi-natedwhosedistributiontocorticalregionsisprovidedinTable1from[13 ],andtherewere71reachingactionsformonkey-1and65formonkey-2,respectively.Thesereach-ingmovementsconsistofthreenaturalsegmentsshowninFigure1. BasedontheanalysisofWessbergetal.[1 ],theinstanta-neousmovementiscorrelatedwiththecurrentandthepastneuraldataupto1second10bins).Therefore,foreachtimeinstant,weformabin-countvectorbyconcatenating10binsofingcounts(whichcorrespondto10-tapdelaylineinalinearer)fromeveryneuron.Hence,ifx j i )rep-resentsthei thbinofneuronj ,wherei 1, ... ,12000} ,a bin-countvectorattimeinstancei isrepresentedbyx i = [ x 1 i ), x 1 i 1), ... x 1 i 9), x 2 i ), ... x n i 9)] T ,wheren is thenumberofneurons.Sinceweareinterestedindetermin-ingrepeatedspatiotemporalingpatternsduringthereach-ingmovements,onlythebincountsfromtimeinstanceswheretheprimatesarmismovingareconsidered.Thereisapossibilitythatintheselectedtrainingsetsomeneuronsneverre.Therowscorrespondingtotheseneuronsmustberemovedfromthebin-countmatrix,sincetheytendtocause 1 AllexperimentalproceduresconformedtotheNationalAcademyPressGuidefortheCareandUseofLaboratoryAnimalsandwereapprovedbytheDukeUniversityAnimalCareandUseCommittee. 010203040506070 30 20 10 0 10 20 30 40 50 60 70 x y z TimemsPosition(mm)ResttofoodFoodtomouthMouthtorest Figure1:Segmentationofthereachingtrajectories:reachfromresttofood,reachfromfoodtomouth,andreachfrommouthtorestpositions(takenfrom[7 ]). instabilityintheNMFalgorithm.Inaddition,topreventtheerrorcriterionfromfocusingtoomuchonneuronsthatsim-plyrefrequently(althoughthetemporalstructureoftheiractivitymightnotbesignicantforthetask),thebincountsineachrowi.e.,foreachneuron)ofthedatamatrixarenor-malizedtohavetheunitlengthinitstwonorms.Ingeneral,if n neuronsareconsideredforatotalofm timeinstances,thedatamatrixX hasdimension(10n m .Sincetheentriesofthedatamatrixarebincounts,theyareguaranteedtobenonnegative.Accountingfor71or65movements,therearem = 2143timeinstancesformonkey-1andm = 2521formonkey-2.3.2.AnalysisoffactorizationprocessIntheapplicationofNMFtoagivenneuralingmatrix,thereareseveralimportantissuesthatmustbeaddressed:theselectionofthenumberofbases,theuniquenessoftheNMFsolution,andunderstandinghowNMFcanndlocalstruc-turesofneuraringactivity.Theproblemofthechoiceofthenumberofbasescanbeaddressedintheframeworkofmodelselection.Anum-berofmodelselectiontechniquese.g.,thecross-validation)canbeutilizedforndingtheoptimalnumberofbases.Inthispaper,wechoosetoadoptaselectioncriterionthathasbeenrecentlydevelopedforclustering.ThecriterioniscalledtheindexI ,whichhasbeenusedtoindicatethecluster

PAGE 4

3116EURASIPJournalonAppliedSignalProcessing validity[15 ].Thisindexhasshownconsistentperformanceofselectingthetruenumberofclustersforvariousexperi-mentalsettings.TheindexI iscomposedofthreefactorsasI r = 1 r E 1 E r D r p ,(4whereE r istheapproximationerror(Frobeniusnormforr bases,andD r isthemaximumEuclideandistancebetweenbasessuchthatD r = r max i j = 1 w i w j (5) Theoptimalr istheonethatmaximizesI r ).Wewillutilizethisindextodeterminetheoptimalr forNMFwithp = 1. DonohoandStoddenhaveshownthatauniquesolu-tionofNMFispossibleundercertainconditions[16 ].TheyhaveshownthroughageometricalinterpretationofNMFthatifthedataarenotstrictlypositive,therecanbeonlyonesetofnonnegativebaseswhichspansthedatainthepositiveorthant.Withanarticulatedsetofimagesobeyingthreerulesagenerativemodel,linearindependenceofgen-erators,andfactorialsampling),theyshowedNMFidenti-thegeneratorsorpartsofimages.Ifweconsiderourneuronalbin-countmatrix,eachrowcontainsmanyzeroen-trieszerobincounts)evenafterremovingnonringneu-ronssincemostneuronsdonotrecontinuouslyonceinev-ery100-millisecondwindowduringtheentiretrainingset.Therefore,ourneuronaldataarenotstrictlypositive.Thisimpliesthattheexistenceofauniquesetofnonnegativebasesfortheneuronalbin-countmatrixiswarranted.Theques-tionstillremainsiftheNMFbasisvectorscanthegen-erativeringpatternsfortheneuralpopulationbymeetingthethreeconditionsmentionedabove.Here,wediscusstheneuronalbin-countdatawithrespecttotheseconditions.Asstatedpreviously,wehavedemonstratedthroughsen-sitivityanalysisthatthespecicneuronalsubsetsfromthePP,PMd,andM1regionsweresequentiallyinvolvedinde-rivingtheoutputofthepredictivemodelsduringreachingmovements[7 ].Hence,thebin-countdataforthereachingmovementwillcontainincreasingingactivityofthespe-cicneuronalsubsetonlocalpartitionsofthetrajectory.Duetobinning,itispossiblethatmorethanoneingpatternisassociatedwithasingledatasample.Thisanalysisleadstoagenerativemodelforthebinneddatainwhichdatasam-plesaregeneratedbylinearcombinationofthespecicringpatternswithnonnegativecoe cients.Also,theseingpat-ternswillbelinearlyindependentsincetheneuronalsubsetineachingpatternstendstomodulateingratesonlyforthelocalpartoftrajectory.Thethirdconditionoffactorialsamplingcanbeapproximatelysatisedbytherepetitionofmovementsinwhichthevariabilityofaparticularingpat-ternisobservedduringtheentiredataset.However,amorerigorousanalysisisnecessarytosupporttheargumentthatthesetofringpatternsiscompleteinfactorialterms.There-fore,weexpectthattheNMFsolutionsmaybeslightlyvari-ableretingtheambiguityinthecompletenessoffactorialsampling.Thismightbeovercomebycollectingmoredataforreachingmovements,andwillbepursuedinfuturestud-ies. 3.3.CasestudiesTheNMFalgorithmisappliedtothedescribedneuronaldatamatrixpreparedusingtentaps,n = 91neuronsformonkey-1aftereliminatingtheneuronsthatdonotethroughtheentiretrainingset)andn = 52neuronsformonkey-2.TheNMFalgorithmwith100independentrunsresultsin r = 5basesforbothmonkey-1andmonkey-2datasetsforwhichtheindexI ismaximized.Themeansandthestan-darddeviationsofthenormalizedcost(Frobeniusnormoferrorbetweenapproximationandthegivendatamatrixdi-videdbytheFrobeniusnormofthedataonly)for100runsare0. 8399 0 001formonkey-1dataand0. 7348 0 002 formonkey-2data.Thisimpliesthatthealgorithmapprox-imatelyconvergestothesamesolutionwithdi erentinitialconditions(althoughnotsu cient). In Figure2,weshowtheresultingbasisvectors(columnsof W forthebincountspresentedinmatrixformwherecolumnsaredi erentneuronsandrowsaredi erentdelays),aswellastheircorrespondingtime-varyingencodingcoe cients(rowsofH superimposedonthereachingtrajectorycoordinatesofthreeconsecutivemovements.Basedontheassumptionthattheneuronalbin-countdataapproximatelysatisfythethreeconditionsfortheidenticationofthegener-ators,theNMFbasisvectorsdeterminethesequenceofspa-tiotemporalringpatternsrepresentingtheingmodula-tionofthespecicneuronalsubsetsduringthecourseofthereachingmovement.Alternatively,wecansaythatNMFdis-coverstheselatentingpatternsofneuralpopulationbyop-timallinearapproximationofthedatawithfewbases[9 ].Forexample,fromthetwobasisvectorseachcorrespondingtotwoprimatesintheleftpanelofFigure2,weobservethatr-ingsoftheneuronsingroup-b arefollowedbyringsoftheneuronsingroup-a (thebrightactivitydenotedbyb occursearlierintimethantheactivitydenotedbya ,sinceincreas-ingvaluesintheverticalaxisofeachbasisindicatesgoingfurtherbackintime).Thus,NMFe ectivelydeterminesandsummarizesthissparseringpatternthatinvolvesagroupofneuronsringsequentially.Theirrelativeaverageactivityisalsoindicatedbytherelativemagnitudesoftheentriesofthisparticularbasis.Usingthesetime-synchronizedneuralactivityandhandtrajectoryrecordings,itisalsopossibletodiscoverrela-tionshipsbetweenringpatternsandcertainaspectsofthemovement.Wecanassesstherepeatabilityofacertainr-ingpatternsummarizedbyabasisvectorbyobservingthetime-varyingactivityofthecorrespondingencodingsignal(thecorrespondingrowofH intime.Anincreaseinthiscoe cientcorrespondstoalargeremphasistothatbasisinreconstructingtheoriginalneuralactivitydata.IntherightpanelofFigure2,weobservethatallbasesareactivatedreg-ularlyintimebytheircorrespondingencodingsignalsatdi erenttimeinstancesandatdi erentamplitudes).Forex-ample,therstbasisformonkey-1isperiodicallyactivatedtothesameamplitude,whereastheactivationamplitudeof

PAGE 5

DeterminationofNeuralFiringPatternsUsingNMF3117 Time501005010050100501005010010 5 1 10 5 1 10 5 1 10 5 1 10 5 1 a b Lag Neuronindex Time(b) 2550255025502550255010 5 1 10 5 1 10 5 1 10 5 1 10 5 1 a b Lag Neuronindex(a) Figure2:(a)Thevebasesformonkey-1top)andmonkey-2bottom).(b)Theircorrespondingencodingsignalsthicksolidline)overlaidonthe3-dimensionalcoordinatesofthereachingtrajectory(dottedlines)forthreeconsecutiverepresentativereachingtasksseparatedbythedashedlines).Notethattheencodingsignalsarescaledtobeinthesameorderofthemagnitudeofthereachingtrajectoryforthevisualpurpose.thethirdbasisvariesineverymovement,whichmightin-dicateachangeintheroleofthecorrespondingneuronalingpatterninexecutingthatparticularmovement.Theperiodicactivationofencodingsalsoindicatestheburstingnatureofthespatiotemporalrepetitivepatterns.Hence,theNMFbasestendtoencodesynchronousandburstingspa-tiotemporalpatternsofneuralringactivity.FromtheNMFdecomposition,weobservecertainasso-ciationsbetweentheactivitiesofneuronsfromdi erentcor-ticalregionsanddi erentsegmentsofthereachingtrajec-tory.Inparticular,ananalysisofthemonkey-1databasedonFigure2indicatesthatneuronsinPPandM1array1repeatsimilarringpatternsduringthereachfromresttofood.Thisassessmentisbasedontheobservationthatbasesthree,four,andve,whichinvolveingactivitiesfromneuronsintheseregions,arerepeatedlyactivatedbytheincreasedamplitudeoftheirrespectiveencodingcoe cients.Similarly,neuronsinM1array2arerepeatedlyactivatedduringthereachtoandfromthemouth(basesoneandtwo.Theseobservationsareconsistentwithourpreviousanalysesthatwereconductedthroughtrainedinput-outputmodelssuchastheWienererandRMLP)[7 ]. Table2comparestheneurons,whichwereobservedtohavethehighestsensitivityfromtrainedmodels,andtheneuronsthathavethelargestmagnitudesineachNMFbasis.Thiscomparisonisbasedonmonkey-1dataset.WecanseethatneuronsfromNMFareasubsetofneuronsobtainedfromthesensitivityanal-ysis.ItisalsoworthstatingthatNMFbasisprovidesmoreinformationthanthemodel-basedsensitivityanalysissinceitdeterminesthesynchronousspatiotemporalpatternswhile

PAGE 6

3118EURASIPJournalonAppliedSignalProcessing Table2:Comparisonofimportantneurons(examinedinthemonkey-1dataset). RegionsPPM1(area1)PMdM1(area2) ThehighsensitiveneuronsthroughRMLP4,5,7,22,26,2938,45None93,94Thelargest-magnitudeneuronsinNMFbases7,2945None93,94 Table3:PerformanceevaluationoftheWienererandthemixtureofmultiplemodelsbasedonNMF. CC(x)CCy)CC(z)MSE(x)MSE(y)MSE(z) Monkey-1 Wienerlter0. 57720. 67120. 75740. 48550. 34680. 2460 NMFmixture0. 71470. 70780. 80760. 27110. 27860. 1627 Monkey-2 Wienerlter0. 37370. 43040. 61920. 30500. 74050. 2882 NMFmixture0. 49740. 50410. 69160. 23540. 54000. 2112 thesensitivityanalysisonlydeterminesindividualimportantneurons.Finally,wewouldliketoreiteratethattheanalysispresentedhereissolelybasedonthedata,whichmeansthatthisanalysisdoesnotneedtotrainaspecicmodeltoinves-tigatetheneuralpopulationorganization.3.4.ModelingimprovementforBMIusingNMFWewilldemonstrateasimpleexampleshowingtheimprovedBMIsperformanceinpredictinghandpositionsbyutilizingNMF.Wewillcomparetheperformanceoftwosystems;theWienererdirectlyappliedtotheoriginalspikecountdataandthemixtureofmultiplelinearersbasedontheNMFbasesandencodings.ThestraightWienerlterisdirectlyappliedtotheneuralingdatatoestimatethethreecoordinatesoftheprimateshandposition.TheWienererhasbeenastandardmodelforBMIs,andmanyotherapproacheshavebeencomparedwithit[19 ].Withninedelays,theinputdimensionalityoftheeris910formonkey-1or510formonkey-2discardinginactivenoing)neuralchannels).Thenweaddabiastoeachinputvectortoestimatethey-intercept.TheweightsoftheerareestimatedbytheWiener-HopfequationasW = R 1 P ,(6whereR isa911 911(or511 511formonkey-2)inputcorrelationmatrix,andP isa911 3or511 3formonkey-2)input-outputcross-correlationmatrix.ThemixtureofmultiplemodelsemploystheNMFen-codingsasmixingcoe cients.AnNMFbasisisusedasawindowfunctionforthecorrespondinglocalmodel.There-fore,eachmodelseesagiveninputvectorthroughadi erentwindowandusesthewindowedinputvectortoproducetheoutput.ThentheNMFencodingsareusedtocombineeachmodelsoutputtoproducethenalestimateofthedesiredhandpositionvector.Thiscanbedescribedinthefollowingequation: d c n = K k = 1 h k n z k n T g k c + b k c ,(7whereh k n isanNMFencodingcoe cientforthek thbasisat n thcolumn(i.e.,timeindex),g k c istheweightvectorofthe k thmodelforthec thcoordinatec [ x y z ]),andb k c is they-interceptofthek thmodelforthec thcoordinate.z k n istheinputvectorwindowedbythek thNMFbasis.Itsi th elementisgivenbyz k i n = x i n w k i (8) Here,x i n isthenormalizedingcountoftheneuroni at timeinstancen ,andw k i isthei thelementofthek thNMFbasis. g k c and b k c canbeestimatedbasedontheMSEcrite-rionbyusingofthestochasticgradientalgorithmsuchasthenormalizedleastmeansquareNLMS).TheweightupdateruleoftheNLMSforeachmodelisthengivenbyg k c n +1= g k c n )+ + z k n 2 h k n e c n z k n ), b k c n +1= b k c n )+ + z k n 2 h k n e c n ), (9) where isthelearningrateand isthenormalizationfactor.e c n istheerrorbetweenthec thcoordinateofthedesiredresponseandthemodeloutput.Intheexperiment,wedividedthedatasamplesinto1771trainingsamplesand372testsamplesformonkey-1datasetand1739and782,respectively,formonkey-2dataset.Theparametersaresetas{ K }={ 0 01,1,5} .Theentiretrainingdatasetispresented60timessu cientenoughfortheweightstoconverge.Theperformanceofthemodelisevaluatedonthetestsetbytwomeasures;thecorrelationcoe cientCC)betweendesiredhandtrajectoryandthemodeloutputtrajectory,andthemeansquarederror(MSE)normalizedbythevarianceofthedesiredresponse.Table3presentstheevaluationoftheperformanceoftwosystemsforbothmonkey-1andmonkey-2datasets.Itshowsasignif-icantimprovementingeneralizationperformancewiththemixtureofmodelsbasedonNMFfactorization.Notethatthegeneralperformanceofmodelsforthemonkey-2datasetisworsethanthatforthemonkey-1

PAGE 7

DeterminationofNeuralFiringPatternsUsingNMF3119 dataset.Thereasonsmaycomefrommanyexperimentalvariables.Oneofthemmaybethenumberofelectrodesandthecorrespondingcorticalareas,aswecanseeinTable1that only32electrodeswereimplantedintwoareasformonkey-2,while64electrodesinfourareasformonkey-1.Toquantifytheperformancedi erencebetweentheWienererandthemixtureofmultiplemodels,wecanap-plyastatisticaltestbasedonthemeansquarederror(MSE)performancemetric[17 ].Bymodelingtheperformancedif-ferenceintermsoftheMSEusingshort-timewindowsasanormalrandomvariable,onecanapplythet-testtoquan-tifysignicance.Thist-testwasappliedtobothmodelingoutputsformonkey-1andmonkey-2with = 0 01or = 0 05.Forbothdatasets,thenullhypothesiswasre-jectedwithbothsignicancelevels,resultinginthep -valuesof0. 0023formonkey-1and0. 0007formonkey-2,respec-tively.Therefore,thestatisticaltestoftheperformancedi er-encedemonstratesthatthemixtureofmultiplemodelsbasedonNMFimprovestheperformancesignicantlycomparedtothestandardWienerer.3.5.DiscussionsTheresultspresentedinthepreviouscasestudyarearepre-sentativeexampleofabroadersetofNMFexperimentsper-formedonthisrecording.Selectionofthenumberoftapsandthenumberofbasesr isdependentontheparticu-larstimulusorbehaviorassociatedwiththeneuraldata.Al-thoughwehaveusedamodelselectionmethodoriginallyde-velopedforclustering,anddidnotprovidefulljusticationthatthisindexissuitabletoNMF,themainmotivationistodemonstratethattheproblemofselectingthenumberofbasescanbeaddressedinthecontextofmodelselection.Thiswillbepursuedinfutureresearch.ThenumberofpatternsthatcanbedistinctlyrepresentedbyNMFislimitedbythenumberofbases.Averysmallnum-berofbaseswillleadtothecombinationofmultiplepatternsintoasinglenonsparsebasisvector.Attheotherextreme,averylargenumberofbaseswillresultinthesplittingofapat-ternintotwoormorebases,whichhavesimilarencodingco-e cientsignalsintime.Inthesesituations,thebasesunderconsiderationcanbecombinedintoonebasis.ItisintriguingthatthemixtureofmodelsbasedonNMFgeneralizesbetterthantheWienererdespitethefactthatthemixturecontainsmuchmoremodelparameters.How-ever,eachmodelinthemixturereceivestheinputsprocessedbythesparsebasisvector.Therefore,eachmodellearnsthemappingbetweenonlyaparticularsubsetofneuronsandhandtrajectories,andthee ectivenumberofparametersforeachmodelismuchlessthanthetotalnumberofinputvari-ables.Moreover,furtheroverttingisavoidedbycombiningtheoutputsoflocalmodelsbythesparseencodingsofNMF.4.CONCLUSIONSNonnegativematrixfactorizationisanovelandrelativelynewtoolforanalyzingthedatastructurewhennonnegativ-ityconstraintsareimposed.InBMIstheneuralinputsareprocessedbygroupingtheingsintobincounts.Sincethebincountsarealwayspositive,wehypothesizedthatNMFwouldbeappropriateforanalyzingtheneuralactivity.Theexperimentalresultsandtheanalysispresentedinthispapershowedthatwecouldrepeatedpatternsinneuronalac-tivitythatoccurredinsynchronywiththereachingbehaviorandwasautomaticallyande cientlyrepresentedinasetofsparsebases.Thesparsenessofthebasesindicatesthatonlyasmallnumberofneuronsexhibitrepeatedingpatternsthatareinuentialinreconstructingtheoriginalneuralac-tivitymatrix.Aspresentedin[10 ],NMFprovideslocalbasesoftheobjects,whileprincipalcomponentanalysis(PCA)providesglobalbases.InourpreliminaryexperimentsofPCAforthesamedata,wehaveobservedthatPCAonlyfoundthemostfrequentlyingneurons,whichmaynotberelatedtothebe-havior.Therefore,NMFcanndlocalrepresentationoftheneuralingdata,andthispropertyofNMFcanbemoreef-fectivethanPCAforBMIswhereingactivitiesofdi erentcorticalareasarecollected.LeeandSeunghaveclaimedintheirpaperthatthesta-tisticalindependenceamongtheencodingsofindependentcomponentanalysisICA)forcesthebasistobeholistic[10 ]. And,iflocalpartsoftheneuralactivityoccurtogetheratthesametime,thecomplicateddependenciesbetweentheen-codingswouldnotbecapturedbytheICAalgorithm.How-ever,wehaveobservedthattheNMFencodingsseemtobeuncorrelatedovertheentiremovement.Hence,ICAwithsomenonnegativeconstraintse.g.,nonnegativeICA[18 ], theICAmodelwithnonnegativebasis[19 ],andnonnegativesparsecoding[20 ])mayyieldinterestingencodingsoftheneuralingactivities.Furtherstudieswillpresentthecom-parisonbetweenNMFandtheseconstrainedICAalgorithmsappliedforBMIs.WhileNMFisfoundtobeausefultoolforanalyzingneuraldatatorepeatableactivitypatterns,therearestillseveralissueswhenusingNMFforneuraldataanalysis.Firstly,themethodonlydetectspatternsofactivity,butitisknownthattheinactivityofaneuroncouldoftenindicatere-sponsetoastimulusorcauseabehavior.AnanalysisbasedonNMFwillfailtoidentifysuchneurons.Next,thenontation-arycharacteristicsofneuralactivitieswouldmakeitdi cult forNMFtoxedspatiotemporalringpatterns.Sincetheneuralensemblefunctiontendstochangeoverneuronalspaceandtimesuchthatdi erentspatio-temporalingpat-ternsmaybeinvolvedforthesamebehavioraloutput,wemayhavetocontinuouslyadaptNMFfactorstotrackthosechanges.ThismotivatesustoconsiderarecursivealgorithmofNMF,whichwillenableustoadaptNMFfactorsonline.Itwillbecoveredinthefuturestudy.InourapplicationofNMF,wedemonstratedthattheNMFlearningalgorithmresultedinsimilarFrobeniousnormoftheerrormatrixfor100runsobtainedwithdif-ferentinitialconditions.However,thisdoesnotnecessarilymeanthattheresultedfactorsaresimilarwithsmallvariance.Therefore,weneedtoquantifythesimilarityoftheNMFre-sultswithdi erentinitializations.Analternativeistoemployothermethodstoobtaintheglobalsolutionsuchasgeneticorsimulatedannealingalgorithms.Thiswillbepresentedinafollow-upreport.

PAGE 8

3120EURASIPJournalonAppliedSignalProcessing ACKNOWLEDGMENTSTheauthorswouldliketothankJohanWessbergforcollect-ingthedatausedinthispaper.ThisworkwassupportedbytheDARPAprojectno.N66001-02-C-8022.REFERENCES [1]J.Wessberg,C.R.Stambaugh,J.D.Kralik,etal.,eal-timepredictionofhandtrajectorybyensemblesofcorticalneuronsinprimates,Nature,vol.408,no.6810,pp.361,2000.[2]D.W.MoranandA.B.Schwartz,otorcorticalactivityduringdrawingmovements:populationrepresentationdur-ingspiraltracing,JournalofNeurophysiology,vol.82,no.5,pp.2693,1999.[3]J.K.Chapin,K.A.Moxon,R.S.Markowitz,andM.A.L.Nicolelis,eal-timecontrolofarobotarmusingsimultane-ouslyrecordedneuronsinthemotorcortex,NatureNeuro-science,vol.2,no.7,pp.664,1999.[4]M.D.Serruya,N.G.Hatsopoulos,L.Paninski,M.R.Fellows,andJ.P.Donoghue,rain-machineinterface:instantneuralcontrolofamovementsignal,Nature,vol.416,no.6877,pp.141,2002.[5]J.C.Sanchez,S.-P.Kim,D.Erdogmus,etal.,nput-outputmappingperformanceoflinearandnonlinearmodelsfores-timatinghandtrajectoriesfromcorticalneuronalingpat-terns,inProc.12thIEEEInternationalWorkshoponNeu-ralNetworksforSignalProcessing,pp.139,Martigny,Switzerland,September2002.[6]S.Darmanjian,S.-P.Kim,M.C.Nechyba,etal.,imodalbrain-machineinterfaceformotorcontrolofroboticpros-thetic,inProc.IEEE/RSJInternationalConferenceonIntelli-gentRobotsandSystems(IROS03),vol.4,pp.3612,LasVegas,Nev,USA,October2003.[7]J.C.Sanchez,D.Erdogmus,Y.N.Rao,etal.,nterpret-ingneuralactivitythroughlinearandnonlinearmodelsforbrainmachineinterfaces,inProc.25thAnnualInternationalConferenceoftheIEEEEngineeringinMedicineandBiologySociety,vol.3,pp.2160,Cancun,Mexico,September2003. [8]J.M.Carmena,M.A.Lebedev,R.E.Crist,etal.,earningtocontrolabrain-machineinterfaceforreachingandgraspingbyprimates,PLoSBiology,vol.1,no.2,pp.1,2003.[9]D.D.LeeandH.S.SeungAlgorithmsfornon-negativema-trixfactorization,inAdvancesinNeuralInformationProcess-ingSystems13,T.K.Leen,T.G.Dietterich,andV.Tresp,Eds.,pp.556,MITPress,Cambridge,Mass,USA,2001.[10]D.D.LeeandH.S.Seung,Learningthepartsofobjectsbynon-negativematrixfactorization,Nature,vol.401,no.6755,pp.788,1999.[11]D.Guillamet,M.Bressan,andJ.Vitri` a,Aweightednon-negativematrixfactorizationforlocalrepresentations,inProc.IEEEComputerSocietyConferenceonComputerVisionandPatternRecognition(CVPR01),vol.1,pp.942,Kauai,Hawaii,USA,December2001.[12]A.vellaandM.C.Tresch,odularityinthemotorsys-tem:decompositionofmusclepatternsascombinationsoftime-varyingsynergies,inAdvancesinNeuralInformationProcessingSystems14,T.G.Dietterich,S.Becker,andZ.Ghahramani,Eds.,pp.629,MITPress,Cambridge,Mass,USA,2002.[13]J.C.Sanchez,Fromcorticalneuralspiketrainstobehav-ior:modelingandanalysis,Ph.D.dissertation,DepartmentofBiomedicalEngineering,UniversityofFlorida,Gainesville,Fla,USA,2004.[14]M.A.L.Nicolelis,A.A.Ghazanfar,B.M.Faggin,S.Votaw,andL.M.Oliveira,econstructingtheengram:simulta-neous,multisite,manysingleneuronrecordings,Neuron, vol.18,no.4,pp.529,1997.[15]U.MaulikandS.Bandyopadhyay,Performanceevaluationofsomeclusteringalgorithmsandvalidityindices,IEEETrans.PatternAnal.MachineIntell.,vol.24,no.12,pp.1650,2002. [16]D.DonohoandV.Stodden,Whendoesnon-negativematrixfactorizationgiveacorrectdecompositionintoparts?inAd-vancesinNeuralInformationProcessingSystems16,S.Thrun,L.K.Saul,andB.Sch olkopf,Eds.,pp.1141,MITPress,Cambridge,Mass,USA,2004.[17]S.-P.Kim,J.C.Sanchez,Y.N.Rao,etal.,ACompar-isonofoptimalMIMOlinearandnonlinearmodelsforbrain-machineinterfaces,submittedtoNeuralComputation, 2004. [18]M.Plumbley,Conditionsfornonnegativeindependentcom-ponentanalysis,IEEESignalProcessingLett.,vol.9,no.6,pp.177,2002.[19]L.Parra,C.Spence,P.Sajda,A.Ziehe,andK.-R.M uller,n-mixinghyperspectraldata,inAdvancesinNeuralInforma-tionProcessingSystems12,S.A.Solla,T.K.Leen,andK.-R.M uller,Eds.,pp.942,MITPress,Cambridge,Mass,USA,2000. [20]P.O.Hoyer,on-negativesparsecoding,inProc.12thIEEEInternationalWorkshoponNeuralNetworksforSignalProcessing,pp.557,Martigny,Switzerland,September2002. Sung-PhilKimwasborninSeoul,SouthKorea.HereceivedaB.S.degreefromtheDepartmentofNuclearEngineering,SeoulNationalUniversity,Seoul,SouthKorea,in1994.In1998,heenteredtheDepartmentofElectricalandComputerEngineering,UniversityofFlorida,inpursuitofMas-terofSciencedegree.HejoinedtheCom-putationalNeuroEngineeringLaboratoryasaResearchAssistantin2000.Healsore-ceivedanM.S.degreeinDecember2000fromtheDepartmentofElectricalandComputerEngineering,UniversityofFlorida.From2001,hecontinuedtopursueaPh.D.degreeintheDe-partmentofElectricalandComputerEngineering,UniversityofFloridaundersupervisionofDr.JoseC.Principe.IntheCompu-tationalNeuroEngineeringLaboratory,hehasinvestigatedthede-codingmodelsandtheanalyticalmethodsforbrain-machinein-terfaces. YadunandanaN.RaowasborninMysore,India.HereceivedhisB.E.degreeinelec-tronicsandcommunicationengineeringfromtheUniversityofMysore,India,inAu-gust1997,andtheM.S.andPh.D.degreesinelectricalandcomputerengineeringfromtheUniversityofFlorida,Gainesville,Fla,in2000and2004,respectively.FromMay2000toJanuary2001,heworkedasaDe-signEngineeratGEMedicalSystems,Wis.CurrentlyheisaSeniorEngineeratMotorola,Fla.Hisresearchinterestsincludeadaptivesignalprocessingtheory,algorithmsandanalysis,neuralnetworksforsignalprocessing,andbiomedicalap-plications.

PAGE 9

DeterminationofNeuralFiringPatternsUsingNMF3121 DenizErdogmusreceivedhisB.S.degreesinelectricalengineeringandmathematicsin1997,andhisM.S.degreeinelectricalengineering,withemphasisonsystemsandcontrol,in1999,allfromtheMiddleEastTechnicalUniversity,Turkey.HereceivedhisPh.D.degreeinelectricalengineeringfromtheUniversityofFlorida,Gainesville,in2002.Since1999,hehasbeenwiththeComputationalNeuroEngineeringLabora-tory,UniversityofFlorida,workingwithJosePrincipe.Hiscurrentresearchinterestsincludeinformation-theoreticaspectsofadaptivesignalprocessingandmachinelearning,aswellastheirapplicationstoproblemsincommunications,biomedicalsignalprocessing,andcontrols.HeistherecipientoftheIEEESPS2003YoungAuthorAward,andisaMemberofIEEE,TauBetaPi,andEtaKappaNu. JustinC.SanchezreceivedaB.S.degreewithhighesthonorsinengineeringsci-encealongwithaminorinbiomechan-icsfromtheUniversityofFloridain2000.From1998to2000,hespentthreeyearsasaResearchAssistantintheDepartmentofAnesthesiology,UniversityofFlorida.In2000,hejoinedtheDepartmentofBiomed-icalEngineeringandComputationalNeu-roEngineeringLaboratory,theUniversityofFlorida.Inthespringof2004,hecompletedbothhisM.E.andPh.D.degreesinbiomedicalsignalprocessingworkingonthede-velopmentofmodelingandanalysistoolsforbrain-machinein-terfaces.HeiscurrentlyaResearchAssistantProfessorintheDe-partmentofPediatrics,DivisionofNeurology,theUniversityofFlorida.Hisneuralengineeringelectrophysiologylaboratoryiscur-rentlydevelopingneuroprostheticsforuseintheresearchandclin-icalsettings. MiguelA.L.NicoleliswasborninSaoPaulo,Brazil,in1961.HereceivedhisM.D.andPh.D.degreesfromtheUniversityofSaoPaulo,Brazil,in1984and1988,re-spectively.AfterpostdoctoralworkatHah-nemannUniversity,Philadelphia,hejoinedDukeUniversity,wherehenowcodirectstheCenterforNeuroengineeringandisaProfessorofneurobiology,biomedicalen-gineering,andpsychologicalandbrainsci-ences.Hislaboratoryisinterestedinunderstandingthegeneralcomputationalprinciplesunderlyingthedynamicinteractionsbe-tweenpopulationsofcorticalandsubcorticalneuronsinvolvedinmotorcontrolandtactileperception. JoseC.PrincipeisaDistinguishedProfes-sorofelectricalandcomputerengineeringandbiomedicalengineeringattheUniver-sityofFloridawhereheteachesadvancedsignalprocessing,machinelearning,andar-ticialneuralnetworks(ANNs)modeling.HeisaBellSouthProfessorandtheFounderandDirectoroftheUniversityofFloridaComputationalNeuroEngineeringLabora-tory(CNEL).Hisprimaryareaofinterestisprocessingoftime-varyingsignalswithadaptiveneuralmodels.TheCNELhasbeenstudyingsignalandpatternrecognitionprin-ciplesbasedoninformation-theoreticcriteria(entropyandmutualinformation).HeisanIEEEFellow.HeisaMemberoftheAD-COMoftheIEEESignalProcessingSociety,aMemberoftheBoardofGovernorsoftheInternationalNeuralNetworkSociety,andtheEditor-in-ChiefoftheIEEETransactionsonBiomedicalEngineer-ing.HeisaMemberoftheAdvisoryBoardoftheUniversityofFloridaBrainInstitute.Hehasmorethan90publicationsinref-ereedjournals,10bookchapters,and200conferencepapers.Hedirected35Ph.D.dissertationsand45Mastertheses.HerecentlywroteaninteractiveelectronicbookentitledNeuralandAdaptiveSystems:FundamentalsThroughSimulationpublishedbyJohnWi-leyandSons.