ManyquestionsinneurosciencedependontheanalysisofneuronalspikingactivityrecordedundervariousbehavioralconditionsForthisreasondataacquiredsimultaneouslyfrommultipleneuronsareinvaluableforeluc ID: 176433
Download Pdf The PPT/PDF document "CommunicatedbyManeeshSahaniUnsupervisedS..." is the property of its rightful owner. Permission is granted to download and print the materials on this web site for personal, non-commercial use only, and to display it on your personal computer provided you do not modify the materials and that you retain all copyright notices contained in the materials. By downloading content from our website, you accept the terms of this agreement.
CommunicatedbyManeeshSahaniUnsupervisedSpikeDetectionandSortingwithWaveletsandSuperparamagneticClusteringR.QuianQuirogarodri@vis.caltech.eduZ.NadasdyDivisionofBiology,CaliforniaInstituteofTechnology,Pasadena,CA91125,U.S.A.Y.Ben-ShaulICNC,HebrewUniversity,Jerusalem,IsraelThisstudyintroducesanewmethodfordetectingandsortingspikesfrommultiunitrecordings.Themethodcombinesthewavelettransform,whichlocalizesdistinctivespikefeatures,withsuperparamagneticclustering,whichallowsautomaticclassiÞcationofthedatawithoutassumptionssuchaslowvarianceorgaussiandistributions.Moreover,animprovedmethodforsettingamplitudethresholdsforspikedetectionisproposed.Wedescribeseveralcriteriaforimplementationthatrenderthealgorithmunsupervisedandfast.Thealgorithmiscomparedtootherconventionalmethodsusingseveralsimulateddatasetswhosecharacteristicscloselyresemblethoseofinvivorecordings.Forthesedatasets,wefoundthattheproposedalgorithmoutperformedconventionalmethods.1Introduction Manyquestionsinneurosciencedependontheanalysisofneuronalspik-ingactivityrecordedundervariousbehavioralconditions.Forthisreason,dataacquiredsimultaneouslyfrommultipleneuronsareinvaluableforelu-cidatingprinciplesofneuralinformationprocessing.Recentadvancesincommerciallyavailableacquisitionsystemsallowrecordingsofuptohun-dredsofchannelssimultaneously,andthereliabilityofthesedatacriti-callydependsonaccuratelyidentifyingtheactivityofindividualneurons.However,developmentsofefÞcientandreliablecomputationalmethodsforclassifyingmultiunitdata,thatis,spikesortingalgorithms,lagbehindthecapabilitiesaffordedbycurrenthardware.Inpractice,supervisedspikesortingofalargenumberofchannelsishighlytime-consumingandnearlyimpossibletoperformduringthecourseofanexperiment.NeuralComputation16,1661Ð16872004MassachusettsInstituteofTechnology 1662R.Quiroga,Z.Nadasdy,andY.Ben-ShaulThebasicalgorithmicstepsofspikeclassiÞcationareasfollows:(1)spikedetection,(2)extractionofdistinctivefeaturesfromthespikeshapes,and(3)clusteringofthespikesbythesefeatures.SpikesortingmethodsaretypicallybasedonclusteringpredeÞnedspikeshapefeaturessuchaspeak-to-peakamplitude,width,orprincipalcomponents(Abeles&Goldstein,1977;Lewicki,1998).Nevertheless,itisimpossibletoknowbeforehandwhichofthesefeaturesisoptimalfordiscriminatingbetweenspikeclassesinagivendataset.InthespeciÞccasewherethespikefeaturesareprojec-tionsontheÞrstfewprincipalcomponents,theplanesontowhichthespikesareprojectedmaximizethevarianceofdatabutdonotnecessarilyprovideanoptimalseparationbetweentheclusters.AsecondcriticalissueisthatevenwhenoptimalfeaturesfromagivendatasetareusedforclassiÞcation,thedistributionofthedataimposesadditionalconstraintsontheclusteringprocedure.Inparticular,violationofnormalityinagivenfeatureÕsdistribu-tioncompromisesmostunsupervisedclusteringalgorithms,andthereforemanualclusteringofthedataisusuallypreferred.However,besidesbeingaverytime-consumingtask,manualclusteringintroduceserrorsduetoboththelimiteddimensionalityoftheclustercuttingspaceandhumanbiases(Harris,Henze,Csicsvari,Hirase,&Buzs«aki,2000).AnalternativeapproachistodeÞnespikeclassesbyasetofmanuallyselectedthresholds(windowdiscriminators)orwithspiketemplates.AlthoughthisiscomputationallyveryefÞcientandcanbeimplementedon-line,itisreliableonlywhenthesignal-to-noiseratioisveryhighanditislimitedtothenumberofchannelsahumanoperatorisabletosupervise.Inthisarticle,weintroduceanewmethodthatimprovesspikesepara-tioninthefeaturespaceandimplementsanovelunsupervisedclusteringalgorithm.Combiningthesetwofeaturesresultsinanovelunsupervisedspikesortingsystem.Thecornerstonesofourmethodarethewavelettrans-form,whichisatime-frequencydecompositionofthesignalwithoptimalresolutioninboththetimeandthefrequencydomains,andsuperparam-agneticclustering(SPC;Blatt,Wiseman,&Domany,1996),arelativelynewclusteringproceduredevelopedinthecontextofstatisticalmechanics.Thecompletealgorithmencompassesthreeprincipalstages:(1)spikedetec-tion,(2)selectionofspikefeatures,and(3)clusteringoftheselectedspikefeatures.IntheÞrststep,spikesaredetectedwithanautomaticamplitudethresh-oldonthehigh-passÞltereddata.Inthesecondstep,asmallsetofwaveletcoefÞcientsfromeachspikeischosenasinputfortheclusteringalgorithm.Finally,theSPCclassiÞesthespikesaccordingtotheselectedsetofwaveletcoefÞcients.Westressthattheentireprocessofdetection,featureextraction,andclusteringisperformedwithoutsupervisionandrelativelyquickly.Inthisstudy,wecomparetheperformanceofthealgorithmwithothermeth-odsusingsimulateddatathatcloselyresemblerealrecordings.Therationaleofusingsimulateddataistoobtainanobjectivemeasureofperformance,sincethesimulationsetstheidentityofthespikes. UnsupervisedSpikeDetectionandSorting16632TheoreticalBackground 2.1WaveletTransform.Thewavelettransform(WT)isatime-frequencyrepresentationofthesignalthathastwomainadvantagesoverconven-tionalmethods:itprovidesanoptimalresolutioninboththetimeandthefrequencydomains,anditeliminatestherequirementofsignalstationar-ity.ItisdeÞnedastheconvolutionbetweenthesignalandthewaveletwherearedilated(contracted),andshiftedversionsofauniquewaveletfunction 2tb wherearethescaleandtranslationparameters,respectively.Equa-tion2.1canbeinverted,thusprovidingthereconstructionofTheWTmapsthesignalthatisrepresentedbyoneindependentvariableontoafunctionoftwoindependentvariables.Thisprocedureisredun-dantandinefÞcientforalgorithmicimplementations;therefore,theWTisusuallydeÞnedatdiscretescalesanddiscretetimesbychoosingthesetofparameters,withintegers.Contractedver-sionsofthewaveletfunctionmatchthehigh-frequencycomponents,whiledilatedversionsmatchthelow-frequencycomponents.Then,bycorrelatingtheoriginalsignalwithwaveletfunctionsofdifferentsizes,wecanobtaindetailsofthesignalatseveralscales.Thesecorrelationswiththedifferentwaveletfunctionscanbearrangedinahierarchicalschemecalledmultires-olutiondecomposition(Mallat,1989).ThemultiresolutiondecompositionalgorithmseparatesthesignalintodetailsatdifferentscalesandacoarserrepresentationofthesignalnamedÒapproximationÓ(fordetails,seeMallat,1989;Chui,1992;Samar,Swartz,&Raghveer,1995;QuianQuiroga,Sakow-icz,Basar,&Sch¬urmann,2001;QuianQuiroga&Garcia,2003).Inthisstudyweimplementedafour-leveldecompositionusingHaarwavelets,whicharerescaledsquarefunctions.Haarwaveletswerechosenduetotheircompactsupportandorthogonality,whichallowsthediscrimi-nativefeaturesofthespikestobeexpressedwithafewwaveletcoefÞcientsandwithoutaprioriassumptionsonthespikeshapes.2.2SuperparamagneticClustering.Thefollowingisabriefdescriptionofthekeyideasofsuperparamagneticclustering(SPC),whichisbasedonsimulatedinteractionsbetweeneachdatapointanditsK-nearestneighbors(fordetails,seeBlattetal.,1996;Blatt,Wiseman,&Domany,1997).ThemethodisimplementedasaMonteCarloiterationofaPottsmodel.The 1664R.Quiroga,Z.Nadasdy,andY.Ben-ShaulPottsmodelisageneralizationoftheIsingmodelwhereinsteadofhavingspinswithvalues2,therearedifferentstatesperparticle(Binder&Heermann,1988).TheÞrststepistorepresenttheselectedfeaturesofeachspikebyainan-dimensionalphasespace.TheinteractionstrengthbetweenisthendeÞnedas Kexp xixj 2 isanearestneighborof0elsewhereistheaveragenearest-neighborsdistanceandisthenumberofnearestneighbors.Notethatthestrengthofinteractionbetweennearest-neighborspikesfallsoffexponentiallywithincreasingEuclideandistance,whichcorrespondstothesimilarityoftheselectedfeatures(i.e.,similarspikeswillhaveastronginteraction).Inthesecondstep,aninitialrandomstatefrom1toisassignedtoeach.ThenMonteCarloiterationsarerunfordifferenttemperaturesusingtheWolfalgorithm(Wolf,1989;Binder&Heermann,1988).Blattetal.(1997)usedaSwendnsen-Wangalgorithminstead,butitsimplementationandperformancearebothverysimilar.TheadvantageofbothalgorithmsoversimplerapproachessuchastheMetropolisalgorithmistheirenhancedperformanceinthesuperparamagneticregime(seeBinder&Heermann,1988;Blattetal.,1997,fordetails).ThemainideaoftheWolfalgorithmisthatgivenaninitialconÞgurationofstates,apointisrandomlyselectedanditsstatechangedtoanewstate,randomlychosenbetween1and.Theprobabilitythatthenearestneighborsofwillalsochangetheirstateisgivenby whereisthetemperature(seebelow).Notethatonlythosenearestneigh-borsofthatwereinthesamepreviousstatearethecandidatestochangetheirvaluesto.NeighborsthatchangetheirvaluescreateaÒfrontierÓandcannotchangetheirvalueagainduringthesameiteration.PointsthatdonotchangetheirvalueinaÞrstattemptcandosoifrevisitedduringthesameiteration.Thenforeachpointofthefrontier,weapplyequation2.4againtocalculatetheprobabilityofchangingthestatetofortheirre-spectiveneighbors.Thefrontierisupdated,andtheupdateisrepeateduntilthefrontierdoesnotchangeanymore.Atthatstage,westarttheprocedureagainfromanotherpointandrepeatitseveraltimesinordertogetrepresen-tativestatistics.Pointsthatarerelativelyclosetogether(i.e.,correspondingtoagivencluster)willchangetheirstatetogether.ThisobservationcanbequantiÞedbymeasuringthepoint-pointcorrelationanddeÞningtobemembersofthesameclusterif,foragiventhreshold UnsupervisedSpikeDetectionandSorting1665AsinBlattetal.(1996),weused20states,11nearestneighbors,500iterations,and5.Ithasindeedbeenshownthatclusteringresultsmainlydependonthetemperatureandarerobusttosmallchangesinthepreviousparameters(Blattetal.,1996).Letusnowdiscusstheroleofthetemperature.Notefromequation2.4thathightemperaturescorrespondtoalowprobabilityofchangingthestateofneighboringpointstogether,whereaslowtemperaturescorrespondtoahigherprobabilityregardlessofhowweaktheinteractionis.Thishasaphysicalanalogywithaspinglass,inwhichatarelativelyhightemperature,allthespinsareswitchingrandomly,regardlessoftheirinteractions(para-magneticphase).Atalowtemperature,theentirespinglasschangesitsstatetogether(ferromagneticphase).However,atacertainmediumrangeoftemperatures,thesystemreachesaÒsuperparamagneticÓphaseinwhichonlythosespinsthataregroupedtogetherwillchangetheirstatesimulta-neously.Regardingourclusteringproblem,atlowtemperatures,allpointswillchangetheirstatetogetherandwillthereforebeconsideredasasinglecluster;athightemperatures,manypointswillchangetheirstateindepen-dentlyfromoneanother,thuspartitioningthedataintoseveralclusterswithonlyafewpointsineach;andfortemperaturescorrespondingtothesuperparamagneticphase,onlythosepointsthataregroupedtogetherwillchangetheirstatesimultaneously.Figure1Ashowsatwo-dimensional(2D)exampleinwhich24002Dpointsweredistributedinthreedifferentclusters.Notethattheclusterspartiallyoverlap,theyhavealargevariance,and,moreover,thecentersfalloutsidetheclusters.Inparticular,thedistancebetweenarbitrarilychosenpointsofthesameclustercanbemuchlargerthanthedistancebetweenpointsindifferentclusters.Thesefeaturesrendertheuseofconventionalclusteringalgorithmsunreliable.Thedifferentmarkersrepresenttheout-comeafterclusteringwithSPC.Clearly,mostofthepointswerecorrectlyclassiÞed.Infact,only102of2400(4%)datapointswerenotclassiÞedbe-causetheywereneartheboundariesoftheclusters.Figure1Bshowsthenumberofelementsassignedtoeachgivenclusterasafunctionofthetem-perature.Atlowtemperatures,wehaveasingleclusterwithall2400datapointsincluded.Atatemperaturebetween0.04and0.05,thisclusterbreaksdownintothreesubclusterscorrespondingtothesuperparamagnetictransi-tion.TheclassiÞcationshownintheupperplotwasperformedatAtabout08,weobservethetransitiontotheparamagneticphase,wheretheclustersbreakdownintoseveralgroupswithafewmemberseach.Notethatthealgorithmisbasedon-nearestneighborinteractionsandthereforedoesnotassumethatclustersarenonoverlappingorthattheyhavelowvarianceoragaussiandistribution.3DescriptionoftheMethod Figure2summarizesthethreeprincipalstagesofthealgorithm:(1)spikesaredetectedautomaticallyviaamplitudethresholding;(2)thewavelettrans- 1666R.Quiroga,Z.Nadasdy,andY.Ben-Shaul Figure1:Exampleshowingtheperformanceofsuperparamagneticclustering.(A)Thetwo-dimensionaldatapointsusedasinputs.Thedifferentmarkersrep-resenttheoutcomeoftheclusteringalgorithm.(B)Clustersizevs.temperature.Attemperature0.05,thetransitiontothesuperparamagneticphaseoccurs,andthethreeclustersareseparated.formiscalculatedforeachofthespikesandtheoptimalcoefÞcientsforseparatingthespikeclassesareautomaticallyselected;and(3)theselectedwaveletcoefÞcientsthenserveastheinputtotheSPCalgorithm,andcluster-ingisperformedafterautomaticselectionofthetemperaturecorrespondingtothesuperparamagneticphase.(AMatlabimplementationofthealgorithmcanbeobtainedon-linefromwww.vis.caltech.edu/rodri.)3.1SpikeDetection.SpikedetectionwasperformedbyamplitudethresholdingafterbandpassÞlteringthesignal(300Ð6000Hz,fourpole UnsupervisedSpikeDetectionandSorting1667 B)C) Figure2:Overviewoftheautomaticclusteringprocedure.(A)Spikesarede-tectedbysettinganamplitudethreshold.(B)AsetofwaveletcoefÞcientsrep-resentingtherelevantfeaturesofthespikesisselected.(C)TheSPCalgorithmisusedtoclusterthespikesautomatically. 1668R.Quiroga,Z.Nadasdy,andY.Ben-ShaulbutterworthÞlter).Thethreshold()wasautomaticallysettoThr whereisthebandpass-Þlteredsignalandisanestimateofthestan-darddeviationofthebackgroundnoise(Donoho&Johnstone,1994).Notethattakingthestandarddeviationofthesignal(includingthespikes)couldleadtoveryhighthresholdvalues,especiallyincaseswithhighÞringratesandlargespikeamplitudes.Incontrast,byusingtheestimationbasedonthemedian,theinterferenceofthespikesisdiminished(underthereason-ableassumptionthatspikesamounttoasmallfractionofallsamples).Todemonstratethis,wegeneratedasegmentof10secofbackgroundnoisewithunitstandarddeviation,andinsuccessivesimulations,weaddedadistinctspikeclasswithdifferentÞringrates.Figure3showsthatfornoisealone(i.e.,zeroÞringrate),bothestimatesareequal,butastheÞringrateincreases,thestandarddeviationofthesignal(conventionalestimate)givesanincreasinglyerroneousestimateofthenoiselevel,whereastheimprovedestimatefromequation3.1remainsclosetotherealvalue.Foreachdetectedspike,64samples(i.e.,2.5ms)weresavedforfurtheranalysis.Allspikeswerealignedtotheirmaximumatdatapoint20.Inordertoavoidspikemisalignmentsduetolowsampling,spikemaximaweredeterminedfrominterpolatedwaveformsof256samples,usingcubic3.2SelectionofWaveletCoefÞcients.Afterspikesaredetected,theirwavelettransformiscalculated,thusobtaining64waveletcoefÞcientsforeachspike.Weimplementedafour-levelmultiresolutiondecompositionusingHaarwavelets.Asexplainedinsection2.1,eachwaveletcoefÞcientcharacterizesthespikeshapesatdifferentscalesandtimes.ThegoalistoselectafewcoefÞcientsthatbestseparatethedifferentspikeclasses.Clearly,suchcoefÞcientsshouldhaveamultimodaldistribution(unlessthereisonlyonespikeclass).Toperformthisselectionautomatically,theLillieforsmodiÞcationofaKolmogorov-Smirnov(KS)testfornormalitywasused(Press,Teukolsky,Vetterling,&Flannery,1992).Notethatwedonotrelyonanyparticulardistributionofthedata;rather,weareinterestedindeviationfromnormalityasasignofamultimodaldistribution.Givenadataset,thetestcomparesthecumulativedistributionfunctionofthedatathatofagaussiandistributionwiththesamemeanandvarianceDeviationfromnormalityisthenquantiÞedbyInourimplementation,theÞrst10coefÞcientswiththelargestdeviationfromnormalitywereused.TheselectedsetofwaveletcoefÞcientsprovides UnsupervisedSpikeDetectionandSorting1669 Figure3:Estimationofnoiselevelusedfordeterminingtheamplitudethresh-old.NotehowtheconventionalestimationbasedonthestandarddeviationofthesignalincreaseswiththeÞringrate,whereastheimprovedestimationfromequation3.1remainsclosetotherealvalue.Seethetextfordetails.acompressedrepresentationofthespikefeaturesthatservesastheinputtotheclusteringalgorithm.Overlappingspikes(i.e.,spikesfromdifferentneuronsappearingquasi-simultaneously)introduceoutliersinthedistributionofthewaveletcoefÞ-cientsthatcausedeviationsfromnormalityinunimodal(aswellasmulti-modal)distributions,thuscompromisingtheuseoftheKStestasanestima-tionofmultimodality.Inordertominimizethiseffect,foreachcoefÞcientweonlyconsideredvalueswithin3standarddeviations.3.3SPCandLocalizationoftheSuperparamagneticPhase.OncetheselectedsetofwaveletcoefÞcientsischosen,weruntheSPCalgorithmforawiderangeoftemperaturesspanningtheferromagnetic,superparamag- 1670R.Quiroga,Z.Nadasdy,andY.Ben-Shaulnetic,andparamagneticphases.Inordertolocalizethesuperparamagneticphaseautomatically,acriterionbasedontheclustersizesisused.Theideaisthatforboththeparamagneticandferromagneticphases,temperaturein-creasescanonlyleadtothecreationofclusterswithfewmemberseach.In-deed,intheparamagneticphase(i.e.,hightemperature),theclustersbreakdownintoseveralsmallones,andintheferromagneticphase,therearealmostnochangeswhenthetemperatureisincreased.Incontrast,inthesuperparamagneticphase,increasingthetemperaturecreatesnewclusterswithalargenumberofmembers.Inourimplementation,wevariedthetemperaturefrom0to0.2inin-crementsof0.01andlookedforthehighesttemperatureatwhichaclus-tercontainingmorethan60pointsappeared(notbeingpresentatlowertemperatures).Sinceoursimulationswere60seclong,thismeansthatweconsideredclusterscorrespondingtoneuronswithameanÞringrateofatleast1Hz.Thethresholdof1Hzgaveusoptimalresultsforalloursim-ulations,butitshouldbedecreasedifoneconsidersneuronswithlowerÞringrates.Alternatively,onecanconsiderafractionofthetotalnumberofspikes.Ifnoclusterwithaminimumof60pointswasfound,wekepttheminimumtemperaturevalue.Usingthiscriterion,wecanautomaticallyselecttheoptimaltemperatureforclusterassignments,andthereforethewholeclusteringprocedurebecomesunsupervised.4DataSimulation Simulatedsignalswereconstructedusingadatabaseof594differentav-eragespikeshapescompiledfromrecordingsintheneocortexandbasalganglia.Forgeneratingbackgroundnoise,spikesrandomlyselectedfromthedatabaseweresuperimposedatrandomtimesandamplitudes.Thiswasdoneforhalfthetimesofthesamples.Therationalewastomimictheback-groundnoiseofactualrecordingsthatisgeneratedbytheactivityofdistantneurons.Next,wesuperimposedatrainofthreedistinctspikeshapes(alsopreselectedfromthesamedatabaseofspikes)onthenoisesignalatrandomtimes.Theamplitudeofthethreespikeclasseswasnormalizedtohaveapeakvalueof1.Thenoiselevelwasdeterminedfromitsstandarddeviation,whichwasequalto0.05,0.1,0.15,and0.2relativetotheamplitudeofthespikeclasses.Inonecase,sinceclusteringwasrelativelyeasy,wealsocon-siderednoiselevels0.25,0.30,0.35,and0.4.Spiketimesandidentitiesweresavedforsubsequentevaluationoftheclusteringalgorithm.ThedatawereÞrstsimulatedatasamplingrateof96,000Hz,andbyusinginterpolatedwaveformsoftheoriginalspikeshapes,wesimulatedthespiketimestofallcontinuouslybetweensamples(tomachineprecision).Finally,thedataweredownsampledto24,000Hz.Thisprocedurewasintroducedinordertoimitateactualrecordingconditionsinwhichsamplesdonotnecessarilyfallonthesamefeatureswithinaspike(i.e.,thepeakofthesignaldoesnotnecessarilycoincidewithadiscretesample). UnsupervisedSpikeDetectionandSorting1671Inallsimulations,thethreedistinctspikeshadaPoissondistributionofinterspikeintervalswithameanÞringrateof20Hz.A2msrefractoryperiodbetweenspikesofthesameclasswasintroduced.Notethattheback-groundnoisereproducesspikeshapevariabilityinbiologicalsignals(Fee,Mitra,&Kleinfeld,1996;Pouzat,Mazor,&Laurent,2002).Moreover,con-structingnoisefromspikesensuresthatthisnoisesharesasimilarpowerspectrumwiththespikesthemselves(1/fspectrum).Therealisticsimula-tionconditionsappliedhererendertheentireprocedureofspikesortingmorechallengingthan,forexample,assumingawhitenoisedistributionofbackgroundactivity.Furthercomplicationsofrealrecordings(e.g.,overlap-pingspikes,burstingactivity,movingelectrodes)willbeaddressedinthenextsection.Figure4showsoneofthesimulateddatasetswithanoiselevel0.1.Figure4Adisclosesthethreespikeshapesthatwereaddedtotheback-groundnoise,asshowninFigure4B.Figure4CshowsasectionofthedatainFigure4BinÞnertemporalresolution.Notethevarianceinshapeandamplitudebetweenspikesofthesameclass(identiÞedwithamarkerofthesamegraylevel)duetotheadditivebackgroundnoise.Figure5showsan-otherexamplewithnoiselevel0.15,inwhichclassiÞcationisconsiderablymoredifÞcultthanintheÞrstdataset.Here,thethreespikeclassessharethesamepeakamplitudesandverysimilarwidthsandshapes.Thedifferences 50 100 150 200 250 300 1 0.5 0 0.5 1 Samples 1 2 3 4 5 1.5 1 0.5 0 0.5 1 1.5 Time (sec) 0.92 0.94 0.96 0.98 1 1.02 1 0.5 0 0.5 1 1.5Time (sec) Class 2 Figure4:Simulateddatasetusedforspikesorting.(A)Thethreetemplatespikeshapes.(B)Thepreviousspikesembeddedinthebackgroundnoise.(C)ThesamedatawithamagniÞedtimescale.Notethevariabilityofspikesfromthesameclassduetothebackgroundnoise. 1672R.Quiroga,Z.Nadasdy,andY.Ben-Shaul 50 100 150 200 250 300 0.5 0 0.5 1 Samples 1 2 3 4 5 1 0.5 0 0.5 1 1.5 Time (sec) 0.08 0.1 0.12 0.14 0.16 0.5 0 0.5 1 1.5Time (sec) Class 2 2212Figure5:Anothersimulateddataset.(A)Thethreetemplatespikeshapes.(B)Thepreviousspikesembeddedinthebackgroundnoise.(C)ThesamedatawithamagniÞedtimescale.HerethespikeshapesaremoredifÞculttodifferen-tiate.NoteinthelowerplotthatthevariabilityinthespikeshapesmakestheirclassiÞcationdifÞcult.betweenthemarerelativelysmallandtemporallylocalized.Byaddingthebackgroundnoise,itappearstobeverydifÞculttoidentifythethreespikeclasses(seeFigure5C).Aswiththepreviousdataset,thevariabilityofspikesofthesameclassisapparent.Allthedatasetsusedinthisarticleareavailableon-lineatwww.vis.rodri.5Results Themethodwastestedusingfourgenericexamplesof60seclength,eachsimulatedatfourdifferentnoiselevels,asdescribedintheprevioussection.SincetheÞrstexamplewasrelativelyeasytocluster,inthiscasewealsogeneratedfourextratimeserieswithhighernoiselevels.5.1SpikeDetection.Figures4and5showtwoofthesimulateddatasets.ThehorizontallinesdrawninFigures4BandCand5BandCarethethresholdsforspikedetectionusingequation3.1.Table1summarizestheperformanceofthedetectionprocedureforalldatasetsandnoiselev-els.Detectionperformancesforoverlappingspikes(i.e.,spikepairswithin64datapoints)arereportedseparately(valuesinbrackets).Overlapping UnsupervisedSpikeDetectionandSorting1673Table1:NumberofMissesandFalsePositivesfortheDifferentDataSets. ExampleNumber(NoiseLevel)NumberofSpikesMissesFalsePositives Example1[0.05]3514(785)17(193)711[0.10]3522(769)2(177)57[0.15]3477(784)145(215)14[0.20]3474(796)714(275)10Example2[0.05]3410(791)0(174)0[0.10]3520(826)0(191)2[0.15]3411(763)10(173)1[0.20]3526(811)376(256)5Example3[0.05]3383(767)1(210)63[0.10]3448(810)0(191)10[0.15]3472(812)8(203)6[0.20]3414(790)184(219)2Example4[0.05]3364(829)0(182)1[0.10]3462(720)0(152)5[0.15]3440(809)3(186)4[0.20]3493(777)262(228)2 Notes:Noiselevelisrepresentedintermsofitsstandarddeviationrelativetothepeakamplitudeofthespikes.Allspikeclasseshadapeakvalueof1.Valuesinbracketsareforoverlappingspikes.spikeshamperthedetectionperformancebecausetheyaredetectedassin-gleeventswhentheyappeartoocloseintime.Incomparisonwiththeotherexamples,arelativelylargenumberofspikeswerenotdetectedindataset1forthehighestnoiselevels(0.15and0.2).Thisisduetothespikeclasswithoppositepolarity(class2inFigure4).Infact,settingupanadditionalnegativethresholdreducedthenumberofmissesfrom145to5fornoiselevel0.15andfrom714to178for0.2.Inthecaseoftheoverlappingspikes,thisreductionisfrom360to52andfrom989to134,respectively.Inallothercases,thenumberofundetectedspikeswasrelativelylow.WiththeexceptionoftheÞrsttwonoiselevelsinexample1andtheÞrstnoiselevelinexample3,thenumberoffalsepositiveswasverysmall(lessthan1%).Loweringthethresholdvalueinequation3.1(e.g.,3.5wouldindeedreducethenumberofmissesbutalsoincreasethenumberoffalsepositives.Theoptimaltrade-offbetweennumberofmissesandfalsepositivesdependsontheexperimenterÕspreference,butweremarkthattheautomaticthresholdofequation3.1givesanoptimalvaluefordifferentnoiselevels.Inthecaseofexample1(noiselevel0.05and0.1)andexample3(noiselevel0.05),thelargenumberoffalsepositivesisexclusivelyduetodoubledetections.Sincethenoiselevelisverylowinthesecases,thethresholdisalsolow,andconsequently,thesecondpositivepeakoftheclass3spikeshowninFigure4isdetected.Onesolutionwouldbetotakeahigherthresholdvalue(e.g.,4.5),butthiswouldnotbeoptimalforhigh 1674R.Quiroga,Z.Nadasdy,andY.Ben-Shaulnoiselevels.Althoughdoubledetectionsdecreasetheperformanceofthedetectionalgorithm,itdoesnotrepresentaproblemwhenconsideringthewholespikesortingprocedure.Inpractice,thefalsepositivesshowupasanadditionalclusterthatcanbedisregardedlater.Forfurthertestingoftheclusteringalgorithm,thecompletedatasetofsimulatedspikes(withboththedetectedandtheundetectedones)willbeused.5.2FeatureExtraction.Figure6showsthewaveletcoefÞcientsforspikesinthedatasetshowninFigure4AandFigure5B.Forclarity,waveletcoefÞcientsofoverlappingspikesarenotplotted.CoefÞcientscorrespond-ingtoindividualspikesaresuperimposed,eachrepresentinghowcloselythespikewaveformmatchesthewaveletfunctionataparticularscaleandtime.CoefÞcientsareorganizedindetaillevels(D)andalastapproxima-tion(A),whichcorrespondtothedifferentfrequencybandsinwhichspikeshapesaredecomposed.EspeciallyinFigure6A,weobservethatsomeofthecoefÞcientsclusterarounddifferentvaluesforthedifferentspikeclasses,thusbeingwellsuitedforclassiÞcation.MostofthesecoefÞcientsarecho-senbytheKStest,asshownwithblackmarkers.Forcomparison,the10coefÞcientswithmaximumvariancearealsomarked.ItisclearfromthisÞgurethatcoefÞcientsshowingthebestdiscriminationarenotnecessarilytheoneswiththelargestvariance.Inparticular,themaximumvariancecri-terionmissesseveralcoefÞcientsfromthehigh-frequencyscales(Dthatallowagoodseparationbetweenthedifferentspikeshapes.Figure7disclosesthedistributionofthe10bestwaveletcoefÞcientsfromFigure6B(inthiscase,includingcoefÞcientscorrespondingtooverlappingspikes)usingtheKScriterionversusthemaximumvariancecriterion.ThreewaveletcoefÞcientsoutofthetenselectedusingtheKScriterionshowade-viationfromnormalitythatisnotassociatedwithmultimodaldistribution:coefÞcient42,showingaskeweddistribution,andcoefÞcients19and10,which,inadditiontoskeweddistribution,havesigniÞcantkurtosismainlyduetotheoutliersintroducedbytheoverlappingspikes.Intheremainingcases,theKScriterionselectedcoefÞcientswithamultimodaldistribution.Incontrast,withtheexceptionofcoefÞcient20,thevariancecriterionselectscoefÞcientswithauniformdistributionthathampersclassiÞcation.Forthesamedata,inFigure8weshowthebestthree-dimensional(3D)projectionsofthewaveletcoefÞcientsselectedwiththeKScriterion(Fig-ure8A),thevariancecriterion(Figure8B)andprojectionsoftheÞrstthreeprincipalcomponents(Figure8C).Inallcases,theclusteringwasdoneauto-maticallywithSPCandisrepresentedwithdifferentgraylevels.WeobservethatusingtheKScriterion,itispossibletoclearlyidentifythethreeclus-ters.Incontrast,whenchoosingthecoefÞcientswiththelargestvariance,itispossibletoidentifytwooutofthreeclusters,andwhenusingtheÞrstthreeprincipalcomponents,onlyasingleclusterisdetected(thenumberofclassiÞcationerrorsisshowninTable2,example2,noise0.15).Notealsothattheclustershapescanbequiteelongated,thuschallenginganycluster-