基于關(guān)聯(lián)維數(shù)的神經(jīng)元動作電位特征提取與分類研究

詹躍榮 范影樂* 楊文偉 楊 勇 李 軼（杭州電子科技大學自動化學院，杭州 3008）

2（杭州電子科技大學生物醫(yī)學工程與儀器研究所，杭州 310018）

引言

動作電位（spike）是神經(jīng)元所接收、分析和傳遞的信號載體，是神經(jīng)元網(wǎng)絡(luò)活動的基本表現(xiàn)形式。近年來，隨著植入式多電極陣列（multi-electrode arrays，MEA）記錄技術(shù)的發(fā)展，用細胞外神經(jīng)元動作電位進行記錄和分析已成為神經(jīng)生理學的重要研究手段。然而，MEA中單個電極所采集到的信號往往是電極附近多個神經(jīng)元發(fā)放的動作電位和大量神經(jīng)系統(tǒng)內(nèi)外噪聲的疊加。因此，如何有效地對spike信號進行檢測和模式分類，是研究神經(jīng)元動作電位序列（spike train）信息特征的重要工作，并可為后續(xù)神經(jīng)元信息編碼研究奠定基礎(chǔ)。

動作電位分類過程主要包括信號檢測、特征提取、模式分類三部分。針對動作電位的特征提取，國內(nèi)外學者提出了主成分分析法（principal component analysis， PCA）［1］、小 波 分 析（wavelet analysis，WA）［2－3］和基于相空間（phase space）的特征提?。?］等多種分析方法。目前，主要的研究方法為主成分分析法和小波分析的方法。其中，PCA方法通過樣本協(xié)方差矩陣求取主要特征值，但由于spike信號具有一定的非線性時變特性，因此 PCA可能無法表征spike信號的完整信息?；谛〔ǚ治龅姆椒ㄍㄟ^對spike信號進行多層分解，在時頻域上對spike信號特征進行動態(tài)描述。但是，此方法得到的信號信息是由小波系數(shù)整體合并來反映的，而在實際分類中小波系數(shù)是單獨作用的，很多有用的波形信息可能被丟失。在模式分類階段，通過對特征向量的分析處理，利用不同的數(shù)據(jù)特性進行統(tǒng)計聚類。目前，大多數(shù)的神經(jīng)元動作電位分類方法是基于動作電位的波形差異，包括模板匹配［5］，以及基于神經(jīng)元放電特性的動作電位分類方法［6］。由于神經(jīng)元的電生理活動存在相互影響，且反映神經(jīng)編碼特性的放電時間間隔分布較為復(fù)雜，因此使用上述方法具有較大的局限性。

神經(jīng)元電生理活動通常被視為處于一個不穩(wěn)定的混沌狀態(tài)［7］，而描述混沌機制的一種有效手段就是分形理論；關(guān)聯(lián)維數(shù)是分形維數(shù)的一種，它作為混沌系統(tǒng)的重要特征，反映了系統(tǒng)的復(fù)雜程度，可用于刻畫神經(jīng)元動作電位的特征信息。本研究提出應(yīng)用關(guān)聯(lián)維數(shù)，對動作電位的波形特征進行無監(jiān)督自動提取?？紤]到K均值（K-means）方法具有較好的收斂性，并能夠在一定程度上解決聚類疊加的問題，所以采用這種方法進行特征聚類，最終實現(xiàn)動作電位的無監(jiān)督分類。

1 相空間重構(gòu)及其相關(guān)參數(shù)的確定

根據(jù)嵌入理論，單個時間序列可以重構(gòu)系統(tǒng)相空間，且重構(gòu)的系統(tǒng)相空間與原始系統(tǒng)拓撲等價。如果嵌入維數(shù)m≥2d+1（d為動力系統(tǒng)維數(shù)），則這個動力系統(tǒng)的吸引子空間幾何結(jié)構(gòu)就會被完全打開，嵌入相空間就可以把有規(guī)律的軌跡恢復(fù)出來。

設(shè)植入式多電極陣列某通道采集到一個spike序列，即

式中，xi為采樣值，n為采樣點個數(shù)。

選擇合適的嵌入延遲時間τ和嵌入維m，得到一個m維的嵌入相空間，相空間中的向量可表示為

令N=n－（m－1）τ，則重構(gòu)的多維相空間可表示為

估計延遲時間τ可通過式（4）～式（6）來確定，有式中，rj=jσ/2，σ為給定時間序列的標準差，C（m，r，t） 為關(guān)聯(lián)積分。

在嵌入維數(shù)m的求取過程中，參考偽最近鄰方法（FNN）的思想進行計算，有

式中，xi（m）為m維相空間中的第i個向量，f（i，m）為第i個向量的最近鄰點的下標。

當m大于某個m0時，若F（m）不再明顯地發(fā)生變化并接近于1，則此時m0+1為最小嵌入維數(shù)［9］。

2 關(guān)聯(lián)維數(shù)

關(guān)聯(lián)維數(shù)以相關(guān)性定義，當系統(tǒng)狀態(tài)隨時間變化時，狀態(tài)變量前后的相關(guān)性可以用來表征信號是否有確定的規(guī)律及其程度，所以關(guān)聯(lián)維數(shù)反映了系統(tǒng)的幾何自相似特性。在確定嵌入延遲時間τ和嵌入維數(shù)m后，在重構(gòu)相空間中利用關(guān)聯(lián)積分計算關(guān)聯(lián)維數(shù)。

設(shè)Yi是重構(gòu)的相空間中第i個向量，計算其余N－1個向量與Yi的距離，采用最大模表示Euclidean距離，即

rij=d（Yi－Yj）=

定義關(guān)聯(lián)積分，有

3 關(guān)聯(lián)維數(shù)提取神經(jīng)元動作電位特征

式中，N為相空間代表點（狀態(tài)矢量）的數(shù)目，ε為相空間中給定超小球的半徑，Θ（·）為 Heaviside函數(shù)。

當ε充分小時，關(guān)聯(lián)維數(shù)可定義為

在實際計算關(guān)聯(lián)維數(shù)時，利用 ln（c?（ε））～lnε曲線，通過線性擬合可得到D2。用高斯函數(shù)u（rij，ε）代替階躍函數(shù)，以改善 ln（c（ε））～lnε關(guān)系的光滑性。高斯函數(shù)定義為

選取一組仿真的典型spike數(shù)據(jù)，如圖1所示。仿真數(shù)據(jù)來源于加州理工大學Andersen實驗室提供的C_Easy1_noise01數(shù)據(jù)，其噪聲水平為0.10，采樣頻率為24 kHz，每個動作電位的時間長度為8/3 ms。圖2（a）表示 C_Easy1_noise01數(shù)據(jù)中3類不同的spike信號，每個子圖為20個同類spike信號疊加在一起，這些同類spike信號在波形上具有類似性，因此狀態(tài)變量的前后相關(guān)性應(yīng)具有類似性。圖2（b）表示對應(yīng)的3類 spike信號的 ln（c（ε））～lnε曲線，其中同類 spike信號的ln（c（ε））～lnε曲線相似度極高，則通過線性擬合得到的關(guān)聯(lián)維數(shù)也是十分接近的；但不同類型 spike信號的 ln（c（ε））～lnε曲線存在較大差異，有利于區(qū)分不同類型的 spike信號。這也表明，利用關(guān)聯(lián)維數(shù)對動作電位模式分類是可行的。

圖1 含噪spike序列Fig.1 Spike train with noise

圖2 3類不同spike信號及其相應(yīng)的ln（c（ε））～lnε曲線。（a）3類不同spike信號波形；（b）3類不同spike對應(yīng)的 ln（c（ε））～lnε曲線Fig.2 Three different types of spike signals and corresponding curve of ln（c（ε））～ lnε.（a）three different types of spike signals；（b）ln（c（ε））～ lnεcurve of three types

在求取關(guān)聯(lián)維數(shù)時，首先對采集到的 spike信號序列進行相空間重構(gòu)。以圖2中的3類spike信號為例，計算確定嵌入延遲時間τ=2，嵌入維數(shù)m=12。如圖3所示，當嵌入維數(shù)m＞11時，F(xiàn)（m）不再明顯發(fā)生變化并且接近于1，取最小嵌入維數(shù)m=12。

圖3 嵌入維數(shù)選擇Fig.3 The diagram of choosing embedding dimension

圖4為來源于上述3類spike信號的疊加，3類spike信號總數(shù)為2 990個。圖5為對應(yīng)3類 spike信號的關(guān)聯(lián)維數(shù)分布，可看出3類spike信號的關(guān)聯(lián)維數(shù)分布具有明顯的分層現(xiàn)象。本研究模式分類算法采用K均值方法，對計算得到的spike關(guān)聯(lián)維數(shù)實現(xiàn)無監(jiān)督分類。

圖4 3類spike信號疊加Fig.4 The superposition of three different types of spike

動作電位分類算法有下列步驟。

步驟1：采用閾值法，檢測隱藏在原始信號中的spike 信號 spi（i=1，2，3，…，n），閾值 Thr確定為［3］

圖5 3類spike信號關(guān)聯(lián)維數(shù)分布Fig.5 The distribution of correlation dimension of three different types of spike

式中：k為經(jīng)驗值，視噪聲而定；E為原始信號經(jīng)巴特沃斯帶通濾波器（帶寬為300～6 000Hz）后的信號，median為取中值。

步驟2：確定嵌入延遲時間τ和合適的嵌入維數(shù)m。

步驟3：計算每個spike所對應(yīng)的關(guān)聯(lián)維數(shù)。

步驟4：利用不同類spike的關(guān)聯(lián)維數(shù)具有明顯的分層分布現(xiàn)象，通過K均值實現(xiàn)spike信號的無監(jiān)督分類。

4 數(shù)據(jù)說明

本算法的驗證數(shù)據(jù)分為仿真數(shù)據(jù)和實驗記錄數(shù)據(jù)兩部分。

仿真數(shù)據(jù)來源于加州理工大學Andersen實驗室提供的spike含噪仿真數(shù)據(jù)，每段神經(jīng)元spike信號采樣時長為6 000 ms，采樣頻率為24 kHz，每個spike片段設(shè)置為8/3 ms，即64個采樣點。仿真數(shù)據(jù)提供每個spike信號所屬的神經(jīng)元類別，為本算法的準確率性能評價提供了參照標準。

實驗記錄數(shù)據(jù)來源于本實驗室采集到的大鼠初級運動皮層神經(jīng)元電信號。通過手術(shù)將16通道的多電極陣列（MEA）植入大鼠的初級運動皮層，然后在術(shù)后第三天采集大鼠該皮層神經(jīng)元的自發(fā)放電序列；MEA采集到的spike模擬信號經(jīng)前置放大器（PreAmp）初步放大后，進入Plexon多通道數(shù)據(jù)采集處理器（MAP），在 MAP提供的40 kHz采樣頻率下轉(zhuǎn)換為spike數(shù)字信號，并在PC機SortClient信號采集軟件配合下，記錄微電極陣列各通道采集的spike序列數(shù)據(jù)。實驗中 MAP系統(tǒng)的增益設(shè)置為1 000，采樣時長為8 000 ms，每個 spike片段時長為1.4 ms，即54個采樣點。大鼠初級運動皮層神經(jīng)元spike數(shù)據(jù)采集系統(tǒng)如圖6所示。

圖6 大鼠初級運動皮層神經(jīng)元spike數(shù)據(jù)采集系統(tǒng)Fig.6 Spike data acquisition system of primary motor cortex in rats

5 結(jié)果及其分析

5.1 仿真數(shù)據(jù)結(jié)果及分析

為了驗證關(guān)聯(lián)維數(shù)對spike信號分類的普遍有效性，筆者隨機選擇4組仿真spike信號進行實驗，每組仿真數(shù)據(jù)為3類不同數(shù)目的spike信號疊加，分類結(jié)果統(tǒng)計如表1所示。其中，準確個數(shù)定義為實驗中每一類正確分類的 spike信號個數(shù)，錯分個數(shù)定義為錯分到該類但實屬它類的spike信號個數(shù)，準確率定義為每類正確分類之和與待分類spike信號的總數(shù)之百分比，1#、2#、3#分別表示每組數(shù)據(jù)中3類不同的spike信號，M和SD分別定義為每組數(shù)據(jù)3類不同spike信號分類準確率的平均值和標準差。從表1中可知，利用關(guān)聯(lián)維數(shù)對spike信號模式進行分類，準確率達到95％以上，每類只有小部分出現(xiàn)了錯分。如前所述，由于受微電極陣列技術(shù)的限制，植入式多電極陣列中的各個電極通常并不能準確定位到單個神經(jīng)元。采集的信號一般來源于神經(jīng)細胞間質(zhì)，是電極周圍多個神經(jīng)元電活動的綜合反映；信號記錄將不可避免地受到外界刺激擾動和周圍神經(jīng)細胞的相互電磁等外干擾的影響，以及細胞體膜參數(shù)、放電閾值的隨機起伏變化等因素引起的內(nèi)噪聲影響。因此，對于部分數(shù)據(jù)，關(guān)聯(lián)維數(shù)方法也出現(xiàn)了錯分，但從表1可看出關(guān)聯(lián)維數(shù)能有效刻畫不同類spike信號的特征，并可達到較為滿意的分類效果。

表1 4組仿真信號的spike分類結(jié)果統(tǒng)計Tab.1 Clustering result of four groups of simulated spikes

表2給出了分別采用PCA方法、小波分析方法和基于關(guān)聯(lián)維數(shù)方法的分類結(jié)果對比，其中PCA方法是利用matlab統(tǒng)計工具箱里的princomp函數(shù)來提取 spike的特征信息，而小波分析方法是利用matlab小波分析工具箱里的wavedec函數(shù)進行小波分解，進而提取spike的特征信息。表2中的1#、2#、3#分別表示每組數(shù)據(jù)中3類不同的spike信號，M和SD分別定義為每組數(shù)據(jù)3類不同spike信號分類準確率的平均值和標準差，其中類間準確率是每一組數(shù)據(jù)中一類spike聚類正確率。從表2中每組數(shù)據(jù)分類的總準確率可得出，基于關(guān)聯(lián)維數(shù)的分類方法總體優(yōu)于PCA分類方法。因類間準確率反映了每類spike信號的分類效果，根據(jù)表2的類間準確率可得出，除個別spike信號PCA方法的分類效果好于基于關(guān)聯(lián)維數(shù)方法的分類效果外，大部分類間準確率顯示，基于關(guān)聯(lián)維數(shù)分類方法優(yōu)于PCA分類方法。表2展示針對不同的數(shù)據(jù)組，小波分析方法和基于關(guān)聯(lián)維數(shù)的分類方法各有長處，總體分類性能大體相當。但小波分析方法針對不同的數(shù)據(jù)，需要選擇合適的小波基來分析，如果面對大數(shù)據(jù)量時，如何選擇合適的小波基將耗費大量時間，影響分類的效率，因此基于關(guān)聯(lián)維數(shù)的分類方法在可實現(xiàn)性方面要優(yōu)于小波分析方法。

5.2 實驗記錄數(shù)據(jù)結(jié)果及分析

微電極陣列第11通道采集的大鼠初級運動皮層神經(jīng)元動作電位序列數(shù)據(jù)，如圖7所示。經(jīng)閾值檢測和濾波后，spike信號的疊加如圖8所示。在對spike信號序列重構(gòu)相空間時，取嵌入延遲時間τ=2，取嵌入維數(shù) m=16。

表2 仿真數(shù)據(jù)的不同分類方法的結(jié)果Tab.2 Clustering result of simulated data wi th different sorting techniques

圖7 大鼠初級運動皮層神經(jīng)元spike序列Fig.7 Spike train of primary motor cortex in rats

從圖8中可見，雖然采集到的大鼠初級運動皮層神經(jīng)元spike信號的波形比較相似，但仍可以明顯看出存在兩類spike信號。圖9為兩類spike信號在相空間中計算出的關(guān)聯(lián)維數(shù)分布，顯示兩類spike信號的關(guān)聯(lián)維數(shù)具有明顯的分層分布現(xiàn)象。圖10為經(jīng)K均值分類后的spike信號疊加，可看出有一部分spike信號錯分，但大部分spike信號分類是正確的。這說明，基于關(guān)聯(lián)維數(shù)的神經(jīng)元動作電位模式分類方法應(yīng)用于真實數(shù)據(jù)是有效的。

圖8 大鼠初級運動皮層神經(jīng)元spike信號疊加Fig.8 Spikes of rats’primary motor cortex

圖9 spike信號關(guān)聯(lián)維數(shù)分布Fig.9 The distribution of spike’s correlation dimension

由于神經(jīng)元動作電位信號是時變、非平穩(wěn)信號，傳統(tǒng)的時域或頻域分析都無法準確地表達神經(jīng)元動作電位的特征信息。基于關(guān)聯(lián)維數(shù)的特征提取優(yōu)勢在于它刻畫了信號的非線性特性，通過將動作電位的時間序列重構(gòu)到高維的相空間中，使得原本在低維空間中無法區(qū)別的特征信息在高維相空間中可以進行有效區(qū)分，增強了非同類動作電位的可區(qū)分度。真實數(shù)據(jù)的實驗結(jié)果表明，本算法框架具有較好的分類效果，代替人工分類具有一定的可行性。

圖10 分類后spike信號疊加Fig.10 The superposition of classified spike

6 結(jié)語

植入式多電極陣列應(yīng)用中的動作電位模式分類技術(shù)，是神經(jīng)編碼研究的前期基礎(chǔ)，因此有著較好的研究意義和應(yīng)用前景。筆者針對神經(jīng)元動作電位模式分類，提出了一種基于關(guān)聯(lián)維數(shù)的特征提取及分類新算法；給出了關(guān)聯(lián)維數(shù)計算過程中延遲時間和嵌入維數(shù)等關(guān)鍵參數(shù)的選擇依據(jù)，以及算法實現(xiàn)的主要步驟。以仿真實驗數(shù)據(jù)為例，結(jié)合K均值模式聚類方法，實現(xiàn)了spike信號的無監(jiān)督分類。與PCA和小波變換等被廣泛使用的方法進行了比較，結(jié)果證明本算法框架對仿真實驗數(shù)據(jù)模式分類的準確率達到了95％以上，表明本算法具有一定的優(yōu)勢。在實驗室環(huán)境下，進行了大鼠初級運動皮層神經(jīng)元電活動的采集，并將本方法應(yīng)用于真實的實驗數(shù)據(jù)，取得了良好效果，證實了非同類神經(jīng)元動作電位的關(guān)聯(lián)維數(shù)具有較好的可分性，同時也驗證了本算法的有效性。

［1］AdamosDA， Kosmidis EK， TheophilidisG. Performance evaluation of PCA-based spike sorting algorithms［J］.Comput Methods Prog Biol，2008，91（3）：232 －244.

［2］Geng Xinling，Hu Guangshu，Tian Xin.Neural spike sorting using mathematical morphology，multiwavelets transform and hierarchical clustering［J］.Neurocomputing，2010，73：707 －715.

［3］Quiroga RQ，Nadasdy Z，Ben-Shaul Y.Unsupervised spike detection and sorting with wavelets and superparamagnetic clustering［J］.Neural Comput，2004，16：1661 － 1687.

［4］Chibirovaa OK， Aksenovaa TI， Benabida A， et al.Unsupervised spike sorting of extracelluar electrophysiological recording in subthalamic nucleus of Parkinsonian patients［J］.Biosystems，2005，79：159－171.

［5］Vargas-Irwin C，Donoghue JP.Automated spike sorting using density grid contour clustering and subtractive waveform decomposition［J］.Journal of Neuroscience Methods，2007，164（1）：1－18.

［6］Delesclue M，Pouzat C.Efficient spike-sorting of multi-state neurons using inter-spike intervals information［J］.Journal of Neuroscience Methods，2006，150：16－29.

［7］彭建華，劉延柱.腦科學中若干非線性動力學問題［J］.力學進展，2003，33（3）：325－321.

［8］Kim HS，Eykholt R，Salas JD.Nonlinear dynamics，delay times and embedding windows［J］.Physica D（S0167 － 2789），1999，127：48－60.

［9］Cao Liangyue.Practical method of determining the minimum embedding dimension of a scalar time series［J］.Physica D，1997，110：43－50.