獨(dú)立分量分析在神經(jīng)元鋒電位分類中的一種新應(yīng)用

楊彭舉 封洲燕 孫 靜

(浙江大學(xué) 生物醫(yī)學(xué)工程與儀器科學(xué)學(xué)院 生物醫(yī)學(xué)工程教育部重點(diǎn)實(shí)驗(yàn)室，杭州 310027)

獨(dú)立分量分析在神經(jīng)元鋒電位分類中的一種新應(yīng)用

楊彭舉 封洲燕*孫 靜

(浙江大學(xué) 生物醫(yī)學(xué)工程與儀器科學(xué)學(xué)院 生物醫(yī)學(xué)工程教育部重點(diǎn)實(shí)驗(yàn)室，杭州 310027)

爆發(fā)式鋒電位(Burst)是大腦神經(jīng)元?jiǎng)幼麟娢话l(fā)放的一種常見形式，它在增強(qiáng)神經(jīng)信號(hào)傳遞的可靠性以及形成突觸可塑性變化等方面具有重要的作用。在細(xì)胞外記錄的鋒電位信號(hào)中，同源Burst也表現(xiàn)為幅值和波形都明顯變化的一串高頻發(fā)放序列，這給神經(jīng)元序列的正確分類提出了難題。為了解決這個(gè)問題，本研究設(shè)計(jì)了一種四極電極記錄的鋒電位信號(hào)檢測(cè)和分類方法。在閾值法檢出鋒電位的基礎(chǔ)上，首先根據(jù)鋒電位時(shí)間間隔指標(biāo)檢出候選Burst信號(hào)小段，然后利用獨(dú)立分量分析(ICA)的盲源分離特性，區(qū)分每個(gè)小段信號(hào)中所包含的不同來(lái)源的鋒電位，再用于最后的整體信號(hào)的鋒電位聚類。實(shí)驗(yàn)記錄數(shù)據(jù)和仿真數(shù)據(jù)的檢驗(yàn)結(jié)果表明，該方法不僅能夠?qū)?lái)自不同神經(jīng)元的Burst和單發(fā)放鋒電位正確分類;而且，由于ICA應(yīng)用的對(duì)象是短時(shí)間的候選Burst信號(hào)，因此，即使對(duì)于4通道信號(hào)也能夠滿足ICA源信號(hào)數(shù)量小于記錄信號(hào)通道數(shù)的限制條件，同時(shí)，短信號(hào)處理又避免了ICA計(jì)算量大等問題，為Burst的正確檢測(cè)與分類提供了一種新方法。

爆發(fā)式發(fā)放;鋒電位;獨(dú)立分量分析;四極電極;分類

引言

大腦神經(jīng)元的動(dòng)作電位發(fā)放模式中，除了單發(fā)放(single spike)以外，還有一種短促的多個(gè)動(dòng)作電位的爆發(fā)式發(fā)放(complex spike burst)，文中簡(jiǎn)稱為Burst。細(xì)胞內(nèi)記錄表明，Burst是由較大的細(xì)胞膜去極化電位誘發(fā)產(chǎn)生，它表現(xiàn)為疊加在同一個(gè)持續(xù)去極化波上的數(shù)個(gè)動(dòng)作電位的高頻率連續(xù)發(fā)放。在細(xì)胞外記錄的動(dòng)作電位(常稱為鋒電位)信號(hào)中，Burst表現(xiàn)為一串幅值和波形都明顯變化的鋒電位［1－2］，脈沖之間的時(shí)間間隔為 2 ～ 20 ms。大腦皮層和海馬組織等區(qū)域的錐體神經(jīng)元都普遍存在Burst形式的動(dòng)作電位發(fā)放，這種發(fā)放在增強(qiáng)神經(jīng)信號(hào)傳遞的可靠性、實(shí)現(xiàn)選擇性信號(hào)傳輸以及產(chǎn)生突觸可塑性效應(yīng)等方面具有重要的作用［3－4］。但是，Burst鋒電位的變化特性給神經(jīng)元發(fā)放序列的檢測(cè)和分析帶來(lái)了難題，同源鋒電位很容易被誤分為來(lái)自不同神經(jīng)元的信號(hào)。至今為止，鋒電位分類的商品化軟件(如 Plexon公司的“offline sorting”等)和“wave_clus”、“MClust”等開源軟件對(duì)于 Burst鋒電位的分類都需要人工干預(yù)，尚不能實(shí)現(xiàn)自動(dòng)正確分類。

利用細(xì)胞外微電極陣列記錄技術(shù)，可以同時(shí)獲得大量神經(jīng)元的鋒電位信號(hào)［5］，這種包含了來(lái)自多個(gè)神經(jīng)元的鋒電位信號(hào)需要進(jìn)行分類，才能夠提取出各個(gè)不同神經(jīng)元的動(dòng)作電位發(fā)放序列［6－7］。常用鋒電位分類的基本依據(jù)是細(xì)胞發(fā)放的動(dòng)作電位具有“全或無(wú)”特性。也就是，如果記錄電極與細(xì)胞之間的距離保持不變，那么，細(xì)胞外某個(gè)特定測(cè)量點(diǎn)上記錄到的來(lái)自同一個(gè)細(xì)胞的鋒電位的幅值和波形應(yīng)該相近，它們的差異是由噪聲引入的［6］。但是，上述Burst所包含的鋒電位卻破壞了這個(gè)原則，因?yàn)閬?lái)自同一個(gè)神經(jīng)元的Burst鋒電位幅值顯著遞減，波形雖然相似也有很大變化，從而難以正確分類［8］。為了解決這個(gè)問題，研究人員提出了許多處理Burst的鋒電位分類算法，大致可以分成兩種:一是先將鋒電位過度分類，然后再將屬于同一個(gè)神經(jīng)元Burst發(fā)放的那些類合并起來(lái);二是根據(jù)Burst中鋒電位發(fā)放的特性先檢出屬于 Burst的鋒電位，然后再進(jìn)行鋒電位分類。

第一種算法的關(guān)鍵是如何將初始的過度分類結(jié)果中屬于Burst的那些類合并起來(lái)。Burst鋒電位波形的漸變特性使它們?cè)诰垲愄卣鲄?shù)空間中表現(xiàn)為拖長(zhǎng)的條帶狀，根據(jù)這個(gè)特點(diǎn)可以挑選出那些屬于同一個(gè)神經(jīng)元的Burst鋒電位［9］。例如，Snider等人將初始分類獲得的每?jī)蓚€(gè)類的所有聚類點(diǎn)投影到穿過這兩類聚類中心點(diǎn)的直線上［10］，如果兩個(gè)聚類中心之間存在的投影點(diǎn)數(shù)量大于某個(gè)設(shè)定的閾值，那么，中間的這些樣本點(diǎn)就被認(rèn)為是Burst中鋒電位漸變所形成的過渡帶，它們與兩個(gè)聚類點(diǎn)一起應(yīng)該合并成同一類。不過，這種算法需要足夠多的用于聚類的鋒電位數(shù)量，否則，即使兩個(gè)聚類屬于同一個(gè)神經(jīng)元，兩個(gè)聚類中心之間的投影點(diǎn)數(shù)量仍然可能達(dá)不到閾值，從而不能將兩個(gè)聚類聯(lián)系起來(lái)。而且，僅僅根據(jù)波形特征參數(shù)空間中聚類點(diǎn)的形式來(lái)合并，很可能誤將不同神經(jīng)元的鋒電位合并在一起。雖然有人通過計(jì)算合并之后的鋒電位發(fā)放時(shí)間間隔(inter-spike interval，ISI)直方圖，檢查直方圖中是否存在清晰的不應(yīng)期，來(lái)判斷合并的正確性［11］，但是，如果初始過度分類所產(chǎn)生的候選合并組合比較多，那么，ISI方法的計(jì)算量就很大。

第二種算法先檢出所有Burst內(nèi)的鋒電位，再取出每個(gè)Burst中的第1個(gè)鋒電位與其它鋒電位一起進(jìn)行分類。這里的關(guān)鍵是如何識(shí)別Burst。如果僅僅依據(jù)Burst內(nèi)鋒電位之間的ISI小于某個(gè)閾值來(lái)確定 Burst［12－14］，顯然，很容易將時(shí)間間隔較小的不同神經(jīng)元的發(fā)放誤判為Burst。如果根據(jù)前后相鄰的鋒電位之間的ISI、幅值衰減系數(shù)以及衰減因子相關(guān)性等多個(gè)參數(shù)來(lái)確定屬于 Burst的鋒電位，可以提高Burst識(shí)別的正確率［8］，但是，這些參數(shù)的取值范圍較難確定，對(duì)于不同的實(shí)驗(yàn)數(shù)據(jù)經(jīng)常需要重新設(shè)置參數(shù)的閾值。

為了解決上述 Burst鋒電位分類的難題，針對(duì)四極電極(tetrode)記錄的多通道信號(hào)設(shè)計(jì)了一種新算法，首先根據(jù)ISI檢出可能包含多個(gè)神經(jīng)元鋒電位的“候選 Burst”;然后，利用獨(dú)立分量分析(independent component analysis，ICA)的盲源分離特性［15］，對(duì)各個(gè)候選Burst的小段信號(hào)進(jìn)行ICA計(jì)算，將其中包含的不同來(lái)源鋒電位分離開來(lái);最后，再進(jìn)行整個(gè)記錄信號(hào)的鋒電位分類。利用ICA處理短時(shí)間的小段信號(hào)是本算法的創(chuàng)新之處，它既能夠滿足ICA源信號(hào)數(shù)量小于記錄信號(hào)通道數(shù)的限制條件，又避免了ICA計(jì)算量大的問題，為Burst的正確檢測(cè)與分類提供了一種新方法。

1 材料和方法

1.1 鋒電位信號(hào)的檢測(cè)和分類方法

圖1所示的四極電極記錄信號(hào)Ch1～Ch4中含有2個(gè)Burst串，分別由5個(gè)和2個(gè)鋒電位構(gòu)成。用該信號(hào)的分析過程來(lái)說明所設(shè)計(jì)的鋒電位信號(hào)的檢測(cè)和分類方法，主要包括4個(gè)步驟。

圖1 包含爆發(fā)式發(fā)放的鋒電位信號(hào)檢測(cè)和分類方法示意Fig.1 Detecting and clustering algorithm for spike signals with bursts

1.1.1 應(yīng)用閾值法檢測(cè)鋒電位

關(guān)鍵是根據(jù)噪聲大小設(shè)定檢測(cè)閾值。由于常用的根據(jù)信號(hào)標(biāo)準(zhǔn)差設(shè)定的閾值受鋒電位發(fā)放率的影響較大，因此，采用基于中值的閾值設(shè)定法［16］，其閾值的計(jì)算公式為

分別計(jì)算每個(gè)通道鋒電位的閾值并檢出各自的鋒電位之后，再取4個(gè)通道檢測(cè)結(jié)果的并集，這樣，就獲得了包含多個(gè)神經(jīng)元鋒電位發(fā)放的時(shí)間序列。

1.1.2 根據(jù)ISI數(shù)值檢出鋒電位序列中的候選Burst

根據(jù)實(shí)驗(yàn)數(shù)據(jù)的觀察結(jié)果，雖然 Burst中相鄰鋒電位之間的ISI都很小，但其變異性較大，很難設(shè)定某個(gè)固定的ISI閾值用于判定不同實(shí)驗(yàn)數(shù)據(jù)的Burst。因此，采用一種自適應(yīng)式計(jì)算 Burst內(nèi)部鋒電位最大ISI值的方法［14］。這種方法對(duì)原始鋒電位時(shí)間序列的ISI數(shù)據(jù)進(jìn)行2次抽取。首先計(jì)算整個(gè)鋒電位發(fā)放序列的平均ISI，表示為

式中，ISIn是第n個(gè)與第n+1個(gè)鋒電位之間的時(shí)間間隔，N為鋒電位總數(shù)。

然后，選出所有小于平均值的 ISIn數(shù)據(jù)，并組成一個(gè)新的ISI序列 L(n)，該 L(n)中包含了所有Burst內(nèi)部的 ISIn［14］。為了進(jìn)一步減少其中的非Burst的ISIn，最后，對(duì) L(n)再求一次平均值，并再抽取其中小于平均值的 ISIn，這些ISIn的平均值就作為候選Burst內(nèi)部鋒電位檢測(cè)的ISI閾值ML。如果ISIn＜ ML，那么，就將其作為 Burst內(nèi)部的 ISIn，再把其中相鄰的 ISIn合并，就完成了候選 Burst的檢測(cè)。

圖1所示的這段鋒電位信號(hào)經(jīng)過2次ISI抽取之后檢出3個(gè)候選Burst。顯然，由于它們只滿足了鋒電位ISI的特性，其中很可能包含了來(lái)自不同神經(jīng)元的鋒電位，因此，下面必須進(jìn)一步對(duì)候選Burst中的鋒電位進(jìn)行分離。

1.1.3 ICA計(jì)算和鋒電位分離

將候選Burst所對(duì)應(yīng)的各個(gè)4通道記錄信號(hào)小段分別進(jìn)行ICA計(jì)算，從而分離出其中包含的不同神經(jīng)元的鋒電位。如果多通道測(cè)量信號(hào)是多個(gè)源信號(hào)的線性組合，并且，各個(gè)源信號(hào)之間統(tǒng)計(jì)獨(dú)立，那么，ICA方法就能夠從測(cè)量信號(hào)中把源信號(hào)分解出來(lái)。來(lái)自不同神經(jīng)元的鋒電位信號(hào)可以滿足這個(gè)要求［15］。設(shè) x=(x1，x2，…xm)T是 m 維測(cè)量信號(hào)，它由源信號(hào) s=(s1，s2，…sn)T中 n個(gè)未知的獨(dú)立源信號(hào)si線性組合而成，即

式中，A為m×n矩陣。ICA分解的關(guān)鍵就是要找到解混矩陣Wn×m，使線性變換 y=Wx求得的向量 y=(y1，y2，…yn)T中的每個(gè)元素 yi(i=1，2，…，n)相互之間都統(tǒng)計(jì)獨(dú)立。

ICA求解 W 的方法有很多，如 FastICA、Infomax、JADE和 RADICAL等算法。比較而言，RADICAL算法根據(jù)統(tǒng)計(jì)獨(dú)立中最自然的判據(jù)—最小互信息原則來(lái)判斷統(tǒng)計(jì)獨(dú)立［16］，從而不必估計(jì)信號(hào)的概率密度，而且也不易受野值的影響。同時(shí)，該算法采用一種快速熵估計(jì)算法，提高了計(jì)算效率。經(jīng)過測(cè)試，該算法對(duì)噪聲的分離效果也很好。因此，使用RADICAL算法實(shí)現(xiàn)ICA盲源分離，并使用文獻(xiàn)［16］提供的MATLAB程序。

4通道候選Burst經(jīng)過ICA計(jì)算后得到4路分離的信號(hào)，這4路信號(hào)中究竟哪些是鋒電位信號(hào)通道和噪聲通道，以及這些通道的排列順序，都是不確定的，鋒電位的峰值取向也會(huì)發(fā)生變化(如負(fù)峰變成正峰)，但其中鋒電位出現(xiàn)的時(shí)刻與原信號(hào)中一致;因此，通過計(jì)算鋒電位時(shí)間位置上各個(gè)通道波形的幅值大小和光滑程度(即計(jì)算局部極值個(gè)數(shù))，來(lái)判定包含鋒電位的信號(hào)通道和噪聲通道［17］。與噪聲相比，鋒電位的波形大而光滑，局部極值個(gè)數(shù)較少，而噪聲的局部極值個(gè)數(shù)較多。

1.1.4 鋒電位分類

經(jīng)過1.1.2和1.1.3兩個(gè)步驟的處理之后，鋒電位信號(hào)所包含的Burst已確定，取出各個(gè)Burst中的第1個(gè)鋒電位與其它所有單發(fā)放的鋒電位一起，再進(jìn)行4通道鋒電位分類計(jì)算。采用的分類算法與文獻(xiàn)［18］相似，簡(jiǎn)介如下:首先將同一鋒電位在4個(gè)通道上的波形連接起來(lái)，形成一個(gè)復(fù)合波形，計(jì)算該復(fù)合波形的主成分(principalcomponent analysis，PCA)，并取90%的主成分分量作為聚類的特征值。然后，應(yīng)用近鄰傳播算法(affinity propagation，AP)對(duì)這些特征值進(jìn)行聚類［19］。

AP聚類算法的特點(diǎn)是無(wú)需預(yù)先知道類別數(shù)量，且分類較精細(xì)，大批量數(shù)據(jù)分類時(shí)效率很高［20］。它利用迭代法尋找聚類中心，首先毫無(wú)偏向地將所有樣本點(diǎn)都視為聚類中心，再計(jì)算兩兩樣本點(diǎn)之間的相似度s(i，k)(相似度可以是歐氏距離或者互相關(guān)系數(shù)等)，構(gòu)成一個(gè)樣本集的相似度矩陣。然后分別從樣本點(diǎn)k和樣本點(diǎn)i的角度出發(fā)，建立另外兩個(gè)參數(shù)r(i，k)和a(i，k)，用于判斷這兩個(gè)樣本點(diǎn)相互選擇的合適度。其中r(i，k)的含義是:與 k作為其它樣本點(diǎn)的聚類中心相比，k作為i聚類中心的合適度。而a(i，k)的含義是:在 i的所有潛在聚類中心里，i選擇 k作為聚類中心的合適度。r(i，k)和 a(i，k)的值越大，k成為最終聚類中心的可能性就越大。這樣，通過反復(fù)迭代，計(jì)算每個(gè)樣本點(diǎn)i到所有其它樣本點(diǎn) k的 r(i，k)和 a(i，k)的值，最終得到的聚類中心就是使(r(i，k)+a(i，k))取值最大的樣本點(diǎn)k。可見，這種算法直接將各個(gè)鋒電位聚集成不同的類，不需要訓(xùn)練過程。AP聚類算法的實(shí)現(xiàn)使用文獻(xiàn)［19］提供的Matlab程序。

聚類完成之后，再加上沒有參加聚類的各個(gè)Burst的其它鋒電位，就得到了圖1最下方所示的最終鋒電位分類結(jié)果。

1.2 鋒電位信號(hào)實(shí)驗(yàn)數(shù)據(jù)的采集

鋒電位信號(hào)采用美國(guó)NeuroNexus Technologies公司生產(chǎn)的微電極陣列在麻醉大鼠海馬CA1和CA3區(qū)采集。如圖2所示，電極陣列的每根電極桿上有4個(gè)間距只有25 μm的記錄點(diǎn)，它們模仿傳統(tǒng)金屬絲制作的四極電極，構(gòu)成一個(gè)菱形。這4個(gè)記錄點(diǎn)由于距離很近，因此，一般都能夠同時(shí)記錄到電極附近來(lái)自同一個(gè)神經(jīng)元的鋒電位，并且通常能夠記錄多個(gè)神經(jīng)元的信號(hào)。各通道記錄信號(hào)的頻率范圍設(shè)定為500 Hz～5 kHz，采樣頻率為20 kHz。

圖2 微電極陣列上菱形排列的記錄點(diǎn)Fig.2 The structure of multisite recording electrodes

1.3 鋒電位仿真信號(hào)的建立

由于實(shí)驗(yàn)記錄中所包含鋒電位信號(hào)的確切信息無(wú)法獲得，為了檢驗(yàn)本研究所設(shè)計(jì)的 Burst鋒電位檢測(cè)和分類算法的有效性，建立了10組4通道仿真數(shù)據(jù)，共包含來(lái)自5個(gè)不同神經(jīng)元的鋒電位信號(hào)。其中3個(gè)神經(jīng)元既有單發(fā)放也有Burst發(fā)放，而另外2個(gè)神經(jīng)元只有單發(fā)放。單發(fā)放和Burst的鋒電位波形都是利用實(shí)驗(yàn)記錄數(shù)據(jù)平均求得的模板。各個(gè)Burst發(fā)放中包含的鋒電位數(shù)有2、3、4個(gè)不等，它們出現(xiàn)的概率根據(jù)實(shí)驗(yàn)數(shù)據(jù)設(shè)定在0.08～0.64的范圍之內(nèi)，其ISI的范圍則設(shè)定為6～20 ms。

根據(jù)神經(jīng)元鋒電位發(fā)放的γ分布規(guī)律［13］，首先生成長(zhǎng)度為120 s的8路數(shù)據(jù)，其中 3路為3類Burst序列，另5路為5類單發(fā)放序列(其中3類分別為3個(gè)Burst神經(jīng)元的單發(fā)放);然后，通過一個(gè)基于正態(tài)分布的4×8隨機(jī)矩陣，將這8路數(shù)據(jù)線性組合成4通道仿真信號(hào)，并在每個(gè)通道上疊加高斯白噪聲。這樣獲得的4通道仿真數(shù)據(jù)中分別包含了來(lái)自5個(gè)神經(jīng)元的統(tǒng)計(jì)獨(dú)立信號(hào)源以及噪聲。各個(gè)神經(jīng)元的鋒電位在各通道中的信噪比不等，鋒電位峰峰幅值與噪聲標(biāo)準(zhǔn)差之比的平均值為9.6±1.6，與實(shí)驗(yàn)記錄數(shù)據(jù)相當(dāng)。仿真數(shù)據(jù)中存在某個(gè)神經(jīng)元的Burst內(nèi)包含其它神經(jīng)元的單發(fā)放和其它Burst發(fā)放，以及不同鋒電位重疊等復(fù)雜情況，由于其中包含的各個(gè)鋒電位及其所屬神經(jīng)元已知，因此，可以用于驗(yàn)證算法的有效性。

1.4 鋒電位檢測(cè)的評(píng)價(jià)參數(shù)

設(shè)NREF為實(shí)際鋒電位個(gè)數(shù)，NS為檢測(cè)獲得的鋒電位總數(shù)，NCS為正確檢出的鋒電位個(gè)數(shù)，那么，根據(jù)文獻(xiàn)［21］中的定義，鋒電位的誤檢率和漏檢率分別為

2 結(jié)果

2.1 實(shí)驗(yàn)數(shù)據(jù)分析結(jié)果

用所設(shè)計(jì)的算法分析了多組實(shí)驗(yàn)記錄信號(hào)，表1所示是其中一組長(zhǎng)22.6 s的4通道記錄信號(hào)所包含的鋒電位信息以及分析結(jié)果。經(jīng)過仔細(xì)的人工目測(cè)挑選，統(tǒng)計(jì)得到這組數(shù)據(jù)中共包含6個(gè)不同神經(jīng)元發(fā)放的1 132個(gè)鋒電位，其中2個(gè)神經(jīng)元具有Burst發(fā)放(即神經(jīng)元1和2)，另外4個(gè)神經(jīng)元(即神經(jīng)元3～6)只有單發(fā)放鋒電位。該組數(shù)據(jù)的鋒電位總漏檢率為1.0%，誤檢率為2.6%。表中3列數(shù)據(jù)從左至右分別是人工目測(cè)的鋒電位統(tǒng)計(jì)數(shù)值(這里不妨作為標(biāo)準(zhǔn))、本算法檢出并分類正確的鋒電位數(shù)，以及算法的正確率(即前2個(gè)數(shù)值的百分比)。該組數(shù)據(jù)的最終鋒電位分類總體正確率為94.3%。

表1 某組實(shí)驗(yàn)數(shù)據(jù)的鋒電位數(shù)量和檢測(cè)分類結(jié)果Tab.1 Result of the algorithm for a typical set of experimental recording

圖3是這組實(shí)驗(yàn)記錄數(shù)據(jù)的4種分析結(jié)果比較，圖中所示是第1和第2主成分分量 PC1和 PC2組成的特征空間聚類投影圖。其中，圖3(a)是人工目測(cè)挑選的結(jié)果，圖3(b)是所設(shè)計(jì)算法的分類結(jié)果，可見兩者基本相符。兩圖中，神經(jīng)元1的 Burst鋒電位發(fā)放在特征空間中呈現(xiàn)明顯的條帶狀，表明了這些鋒電位波形的漸變過程。盡管它們?cè)诰垲惪臻g中很分散，但本算法仍然能夠正確分類。如果跳過1.1節(jié)所述的第2和第3個(gè)Burst處理步驟，直接將所有測(cè)得的鋒電位一起進(jìn)行分類，那么，如圖3(c)所示，鋒電位數(shù)量較多的神經(jīng)元1的鋒電位就被分成了4類。

如果在應(yīng)用ISI閾值法提取候選Burst后，就將它們作為Burst，而不用ICA法做鋒電位分離，也就是跳過1.1節(jié)算法的第3步驟，直接取出各個(gè)Burst的第1個(gè)鋒電位參與后續(xù)分類;那么，如圖3(d)所示，由于在神經(jīng)元1的Burst期間或者緊隨其后，夾雜著許多其它單發(fā)放鋒電位(它們主要來(lái)自神經(jīng)元3、4和6)，這些鋒電位在 ISI閾值法提取 Burst時(shí)就直接被誤歸入神經(jīng)元1的Burst中。同時(shí)，圖3(d)還表明有許多神經(jīng)元1的鋒電位被誤歸入神經(jīng)元3、4和6的類別中，這是由于神經(jīng)元1鋒電位之前在ISI閾值范圍內(nèi)出現(xiàn)了這些神經(jīng)元的鋒電位。由此可見ICA分離操作的重要性。

其實(shí)，這組實(shí)驗(yàn)數(shù)據(jù)的鋒電位發(fā)放情況比較復(fù)雜，其中包含的大部分Burst內(nèi)都夾雜著其它神經(jīng)元的單發(fā)放鋒電位，甚至還有2個(gè)不同神經(jīng)元的Burst互相交叉的情況。圖4的兩張圖分別顯示了其中的2個(gè)候選 Burst及其 ICA分離結(jié)果，圖的左邊是4通道記錄數(shù)據(jù)，其上方的短豎線表示檢出的鋒電位及其序數(shù)。比較4路信號(hào)可以看出，圖4(a)中的第1～5個(gè)以及第7個(gè)鋒電位是同一個(gè)神經(jīng)元的Burst發(fā)放，而第6個(gè)鋒電位來(lái)自另一個(gè)神經(jīng)元。經(jīng)過ICA分解之后，這段候選 Burst中屬于2個(gè)神經(jīng)元的鋒電位被分離在第1和2個(gè)通道上，另兩個(gè)通道則是噪聲。

圖4(b)的情況更復(fù)雜，它共包含來(lái)自4個(gè)神經(jīng)元的鋒電位。其中第 1、3、4、6個(gè)鋒電位組成的Burst與圖4a的Burst是同一類，只是鋒電位個(gè)數(shù)不同。而第 5、7、8、9、10 個(gè)鋒電位組成的 Burst則來(lái)自另一個(gè)神經(jīng)元。值得注意的是，圖4b的第3個(gè)鋒電位波形其實(shí)是2個(gè)神經(jīng)元的鋒電位的疊加，相當(dāng)于圖4a中第2個(gè)和第6個(gè)鋒電位的重疊。如果不仔細(xì)分析，這種重疊肉眼很難辨別，而 ICA卻能夠正確地將重疊的2個(gè)鋒電位分開，這進(jìn)一步說明了ICA分析的作用。但是，除了噪聲通道以外，4路ICA信號(hào)最多只能區(qū)分3類鋒電位信號(hào)，因此，圖4(b)中的第2個(gè)鋒電位(屬于第4個(gè)神經(jīng)元)就沒有被分出來(lái)，不過，在如此短的時(shí)間段中像圖4(b)這樣復(fù)雜的鋒電位發(fā)放情況較少見。

圖3 某組實(shí)驗(yàn)數(shù)據(jù)的聚類結(jié)果。(a)根據(jù)人工目測(cè)統(tǒng)計(jì)結(jié)果作出的聚類圖;(b)本算法的聚類圖;(c)未進(jìn)行Burst處理直接將所有鋒電位進(jìn)行分類的結(jié)果;(d)根據(jù)ISI檢出Burst但不使用ICA分離處理而直接聚類的結(jié)果Fig.3 Clustering results of a typical experimental recording.(a)human visual inspection;(b)using the complete algorithm;(c)using the algorithm lack of burst processing;(d)using the algorithm lack of ICA separation

筆者統(tǒng)計(jì)了2 630個(gè) Burst信號(hào)段，其中包含4類以上鋒電位的只占1.4%;因此，絕大多數(shù)候選Burst小段信號(hào)中只含有3類及以下鋒電位，滿足ICA源信號(hào)數(shù)量小于記錄信號(hào)通道數(shù)的限制條件。

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

由于噪聲和不同鋒電位重疊等復(fù)雜因素的影響，實(shí)驗(yàn)記錄中的真實(shí)鋒電位發(fā)放信息無(wú)法獲得，從而無(wú)法精確統(tǒng)計(jì)算法的分析結(jié)果，用10組仿真數(shù)據(jù)檢驗(yàn)了本研究所設(shè)計(jì)算法的正確性，其結(jié)果如表2所示。這10組數(shù)據(jù)的鋒電位平均總數(shù)為1 326±22個(gè)，本算法的鋒電位漏檢率為(2.6±0.48)%，誤檢率為(1.8±0.7)%。仿真數(shù)據(jù)的3個(gè)具有Burst發(fā)放的神經(jīng)元的鋒電位中，屬于Burst的鋒電位總數(shù)分別為(221±14)、(134 ±9)和(72 ±5)個(gè)，單發(fā)放鋒電位總數(shù)分別為(23 ±2)、(51 ±4)和(74±2)個(gè)。Burst檢出率定義為檢出的 Burst鋒電位總數(shù)與Burst實(shí)際包含的鋒電位總數(shù)之比。本算法的Burst平均檢出率達(dá)(95.1±1.7)%，包括單發(fā)放的所有鋒電位分類的最終總體正確率為(96.5±1.0)%。

圖5所示的是某組仿真數(shù)據(jù)的聚類結(jié)果在前3個(gè)主成分特征分量(即PC1、PC2和PC3)組成的特征空間中的投影，其中圖5(a)是根據(jù)已知的仿真信息建立的完全正確的理想分類圖，圖5(b)是文中1.1節(jié)所述算法的分類結(jié)果，而圖5(c)則是跳過該算法的第2個(gè)和第3個(gè)Burst處理步驟，直接將所有檢測(cè)的鋒電位一起進(jìn)行分類的結(jié)果。與實(shí)驗(yàn)數(shù)據(jù)的檢測(cè)結(jié)果相似，圖5(b)與圖5(a)的結(jié)果基本吻合，只有少量誤檢和誤分類。仿真數(shù)據(jù)中包含的3類Burst的鋒電位在特征空間中分別呈現(xiàn)出沿3條直線排列的多個(gè)聚類(見圖5(a)和(b))，而圖5(c)卻把這些屬于同一類Burst的鋒電位分成了不同的類，而且還把不屬于同一類的 Burst誤分成了同一類。

圖4 候選Burst及其ICA分離示例。(a)信號(hào)中包含2類鋒電位和1個(gè)Burst;(b)信號(hào)中包含4類鋒電位和2個(gè)BurstFig.4 Illustration of candidate burst and its ICA separation.(a)a signal segment with 3 types of spike including 1 burst;(b)a signal segment with 4 types of spike including 2 bursts

圖5 某組仿真數(shù)據(jù)的聚類結(jié)果。(a)完全正確的理想分類圖;(b)本算法的分類圖;(c)未進(jìn)行Burst處理直接將所有鋒電位進(jìn)行分類的結(jié)果Fig.5 Clustering result for a typical set of synthetic signals.(a)real spike clusters;(b)clusters of the complete algorithm;(c)clusters of the algorithm without burst detection and ICA separation.Different colors represent different spike clusters

表2 仿真數(shù)據(jù)的鋒電位數(shù)量和檢測(cè)分類結(jié)果(n=10)Tab.2 Results of the algorithm for 10 sets of synthetic signals

3 討論

所設(shè)計(jì)的鋒電位檢測(cè)和分類算法較好地處理了Burst鋒電位難以正確分類的問題，其特點(diǎn)是根據(jù)ISI特性檢出候選 Burst信號(hào)小段之后，先利用ICA盲信號(hào)分離技術(shù)確定Burst小段中包含的不同類鋒電位，然后再應(yīng)用主成分分量對(duì)整個(gè)信號(hào)的鋒電位進(jìn)行分類。其中ICA的應(yīng)用具有如下特點(diǎn):

(1)應(yīng)用 ICA處理候選 Burst信號(hào)小段。ICA分離技術(shù)的獨(dú)特優(yōu)勢(shì)是能夠識(shí)別Burst這樣的來(lái)自同一個(gè)信號(hào)源而波形卻明顯變化的非穩(wěn)態(tài)鋒電位，并且能夠分離幾乎完全重疊的鋒電位(圖4)，這些是其它方法很難做到的［6］。而且，Burst期間，同一個(gè)神經(jīng)元高頻率發(fā)放動(dòng)作電位，其鋒電位與其它神經(jīng)元的鋒電位重疊的概率要比動(dòng)作電位單發(fā)放期間高，這進(jìn)一步體現(xiàn)了本研究算法中應(yīng)用ICA技術(shù)的重要性。

(2)四極電極記錄信號(hào)符合ICA的應(yīng)用要求。ICA只適用于多個(gè)通道能夠同時(shí)記錄到來(lái)自同一信號(hào)源的情況，因此，它常用于腦電圖EEG之類頻率較低、傳播距離較遠(yuǎn)的多通道信號(hào)分析［22］。神經(jīng)元鋒電位(即細(xì)胞外記錄的動(dòng)作電位)由于頻率高、衰減快，相距稍遠(yuǎn)的電極記錄點(diǎn)就不能同時(shí)記錄到同一個(gè)神經(jīng)元的鋒電位，即使是多通道記錄信號(hào)，如果每個(gè)通道記錄到的神經(jīng)元都不相同，也無(wú)法使用ICA分離技術(shù)。不過，實(shí)際上僅僅根據(jù)單通道記錄信號(hào)并不能正確判斷鋒電位是否來(lái)自同一個(gè)神經(jīng)元，因此，鋒電位信號(hào)采集中經(jīng)常使用相距很近的四極電極［9］，如 4根金屬絲并在一起制作成的電極，或者使用本研究的微電極陣列等。這類電極記錄的4通道信號(hào)能夠同時(shí)記錄到相同神經(jīng)元的信號(hào)，可以利用ICA技術(shù)。

(3)候選Burst信號(hào)小段的 ICA分析滿足通道數(shù)的限制條件。ICA分離技術(shù)有一個(gè)限制條件，即源信號(hào)的數(shù)量要小于記錄的通道數(shù)，這就意味著除了噪聲以外，四極電極記錄的4通道信號(hào)最多只能識(shí)別3個(gè)神經(jīng)元的鋒電位。而四極電極能夠記錄到的總的神經(jīng)元數(shù)很容易超過這個(gè)數(shù)量，雖然有人設(shè)計(jì)了許多算法來(lái)克服這個(gè)限制條件［23］，再將ICA用于鋒電位分類，但效果不佳。不過，本研究設(shè)計(jì)的算法卻能夠滿足這個(gè)限制條件，因?yàn)橹粚?duì)各個(gè)候選Burst的小段信號(hào)進(jìn)行 ICA計(jì)算，由于時(shí)間很短，這些小段信號(hào)98%以上包含的鋒電位類別不超過3類。實(shí)驗(yàn)數(shù)據(jù)和仿真數(shù)據(jù)較高的總體鋒電位分類正確率也證明了這一點(diǎn)。

(4)只對(duì)一部分信號(hào)小段進(jìn)行ICA分析也大大減少了ICA計(jì)算量。長(zhǎng)時(shí)間記錄信號(hào)分析時(shí)計(jì)算量過大也是阻礙ICA應(yīng)用的一個(gè)大問題。因此，本研究的算法既發(fā)揮了ICA的優(yōu)勢(shì)又克服了其缺點(diǎn)。

另外，有關(guān)Burst內(nèi)部最大鋒電位時(shí)間間隔 ISI的數(shù)值，Harris等設(shè)定為6 ms，并將 ISI為7～20 ms的稱為偽爆發(fā)［2］。筆者根據(jù)4通道實(shí)驗(yàn)記錄信號(hào)中各個(gè)通道上鋒電位的形態(tài)變化，發(fā)現(xiàn) Burst內(nèi)鋒電位的 ISI變化較大，因此，本研究檢測(cè)候選 Burst時(shí)采用了2次抽取ISI的自適應(yīng)方法來(lái)確定閾值［14］，而不是設(shè)定某個(gè)固定的 ISI閾值。共統(tǒng)計(jì)了25組實(shí)驗(yàn)數(shù)據(jù)，計(jì)算得到兩次“抽取”之后所得ISI的平均值為(17.2 ±10.1)ms。該數(shù)值遠(yuǎn)比文獻(xiàn)中常用的固定 ISI(如6 ms)要大，因此，不會(huì)漏檢Burst。而且，將候選 Burst做進(jìn)一步處理后才獲得真正的Burst，即使候選 Burst中包含了多余的較大ISI的非Burst信號(hào)也無(wú)妨。如果需要檢測(cè)滿足某個(gè)固定ISI閾值的 Burst，那么，可以在檢測(cè)候選 Burst時(shí)直接使用固定閾值，或者在分類結(jié)束之后的鋒電位序列中，再用ISI的固定閾值檢出Burst。

4 結(jié)論

本研究利用ICA盲源分離技術(shù)設(shè)計(jì)的方法不僅能夠?qū)儆贐urst的鋒電位和單發(fā)放鋒電位正確分類，而且還能夠分離重疊鋒電位，同時(shí)又滿足了ICA源信號(hào)數(shù)量的限制條件，并避免了ICA計(jì)算量大等問題，為Burst的正確檢測(cè)與分類提供了一種新方法。

［1］ Roberts MT，Bender KJ，Trussell LO.Fidelity of complex spike-mediated synaptic transmission between inhibitory interneurons［J］.J Neurosci，2008，28(38):9440 － 9450.

［2］ Harris KD，Hirase H，Leinekugel X，et al.Temporal interaction between single spikes and complex spike bursts in hippocampal pyramidal cells［J］.Neuron，2001，32(1):141 － 149.

［3］ Izhikevich EM，Desai NS，Walcott EC，et al.Bursts as a unit of neural information:selective communication via resonance［J］.Trends Neurosci，2003，26(3):161 －167.

［4］ Krahe R，Gabbiani F.Burst firing in sensory systems［J］.Nat Rev Neurosci，2004，5(1):13 － 23.

［5］ 封洲燕，光磊，鄭曉靜，等.應(yīng)用線性硅電極陣列檢測(cè)海馬場(chǎng)電位和單細(xì)胞動(dòng)作電位［J］.生物化學(xué)與生物物理進(jìn)展，2007，34(4):401－407.

［6］ Lewicki MS.A review of methods for spike sorting:the detection and classification of neural action potentials［J］.Network，1998，9(4):R53－R78.

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

［8］ Kaneko H，Tamura H，Suzuki SS.Tracking spike-amplitude changes to improve the quality of multineuronal data analysis［J］.IEEE Trans Biomed Eng，2007，54(2):262 －272.

［9］ Gray CM，Maldonado PE，Wilson M，et al.Tetrodes markedly improve the reliability and yield of multiple single-unit isolation from multi-unit recordings in cat striate cortex［J］.J Neurosci Methods，1995，63(1－2):43－54.

［10］ Snider RK，Bonds AB.Classification of non-stationary neural signals［J］.J Neurosci Methods，1998，84(1 －2):155 －166.

［11］ Fee MS，Mitra PP，Kleinfeld D.Automatic sorting of multiple unit neuronal signals in the presence of anisotropic and non-Gaussian variability［J］.J Neurosci Methods，1996，69(2):175－188.

［12］ 陳琳，林云生，陳傳平，等.基于最大鋒電位間隔的爆發(fā)檢測(cè)自適應(yīng)算法［J］.生物物理學(xué)報(bào)，2006，22(4):318－324.

［13］ Gourevitch B，Eggermont JJ.A nonparametric approach for detection of bursts in spike trains［J］.J Neurosci Methods，2007，160(2):349－358.

［14］ Chen Lin，Deng Yong，Luo Weihua，et al.Detection of bursts in neuronal spike trains by the mean inter-spike interval method［J］.Progress in Natural Science，2009，19(2):229 －235.

［15］ Brown GD，Yamada S，Sejnowski TJ.Independent component analysis at the neural cocktail party［J］. Trends Neurosci，2001，24(1):54－63.

［16］ Learned-Miller EG，F(xiàn)isher JW.ICA using spacings estimates of entropy［J］.Journal of Machine Learning Research，2004，4(7－8):1271－1295.

［17］ Greene BR，F(xiàn)aul S，Marnane WP，et al.A comparison of quantitative EEG features for neonatal seizure detection［J］.Clin Neurophysiol，2008，119(6):1248 －1261.

［18］ 王靜，封洲燕.多通道神經(jīng)元鋒電位檢測(cè)和分類的新方法［J］.生物化學(xué)與生物物理進(jìn)展，2009，36(5):641－647.

［19］ Frey BJ，Dueck D.Clustering by passing messages between data points［J］.Science，2007，315(5814):972 － 976.

［20］ Kelly. Affinity Program SlashesComputing Times［OL］.http://www.scientificblogging.com/cash/affinity_program_slashes_computing_times.2007－10－25/2011－03－10.

［21］ Maccione A，Gandolfo M，Massobrio P，et al.A novel algorithm for precise identification of spikes in extracellularly recorded neuronal signals［J］.J Neurosci Methods，2009，177(1):241－249.

［22］ Jervis BW，Belal S，Cassar T，et al.Waveform analysis of nonoscillatory independent components in single-trial auditory eventrelated activity in healthy subjects and Alzheimer＇s disease patients［J］.Curr Alzheimer Res.2010，7(4):334 － 347.

［23］ Takahashi S，Anzai Y，Sakurai Y.A new approach to spike sorting for multi-neuronal activities recorded with a tetrode-how ICA can be practical［J］.Neurosci Res，2003，46(3):265 －272.

Novel Application of Independent Component Analysis on Neuronal Spike Sorting

YANG Peng-Ju FENG Zhou-Yan*SUN Jing
(College of Biomedical Engineering and Instrumentation Science，Key Lab for Biomedical Engineering of Education Ministry，Zhejiang University，Hangzhou 310027，China)

Complex spike burst(i.e.，Burst)is one of the common modes of neuronal action potential firing in brain.Evidences have showed that burst firings play important roles in increasing reliability of neuronal signal transmission，in generating synaptic plasticity and so on.In extracellular recordings，a burst appears as a train of high-frequency firing spikes with obvious changes both in amplitude and in waveform of spikes.The nonsteady change feature of burst spikes is a challenge to the accurate analysis of neuronal firing sequences.This paper presented a spike sorting algorithm for tetrode recording signals in order to deal with the burst spikes.In this method，based on the spike signals collected by a threshold detecting method，short candidate burst segments were selected firstly according to inter-spike intervals.Then the independent component analysis(ICA)was used to separate spikes from different sources in each of the short signal segments.Finally spike clustering for the whole signals was fulfilled based on the results of ICA.The results obtained from both experimental recordings and synthetic data showed that the algorithm was able to sort the burst spikes and single spikes from different sources accurately.In addition，because the ICA is applied on the very short signal segments of candidate bursts，even for the only four channel signals the algorithm can also meet the source signal number limit of ICA technique.Further more，the short signal processing decreases the ICA computing time.Therefore，the algorithm provides a new method for accurate burst spike detecting and sorting.

complex spike burst;spike;independent component analysis;tetrode;sorting

R338，TP391

A

0258-8021(2011)04-0500-09

10.3969/j.issn.0258-8021.2011.04.004

2011-03-20，錄用日期:2011-05-22

國(guó)家自然科學(xué)基金(30770548，30970753).

* 通信作者。 E-mail:hnfzy@yahoo.com.cn