楊彭舉 封洲燕 孫 靜

(浙江大學 生物醫(yī)學工程與儀器科學學院 生物醫(yī)學工程教育部重點實驗室，杭州 310027)

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

楊彭舉 封洲燕*孫 靜

(浙江大學 生物醫(yī)學工程與儀器科學學院 生物醫(yī)學工程教育部重點實驗室，杭州 310027)

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

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

引言

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

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

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

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

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

1 材料和方法

1.1 鋒電位信號的檢測和分類方法

圖1所示的四極電極記錄信號Ch1～Ch4中含有2個Burst串，分別由5個和2個鋒電位構成。用該信號的分析過程來說明所設計的鋒電位信號的檢測和分類方法，主要包括4個步驟。

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

1.1.1 應用閾值法檢測鋒電位

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

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

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

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

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

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

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

1.1.3 ICA計算和鋒電位分離

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

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

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

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

1.1.4 鋒電位分類

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

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

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

1.2 鋒電位信號實驗數(shù)據(jù)的采集

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

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

1.3 鋒電位仿真信號的建立

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

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

1.4 鋒電位檢測的評價參數(shù)

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

2 結果

2.1 實驗數(shù)據(jù)分析結果

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

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

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

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

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

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

圖3 某組實驗數(shù)據(jù)的聚類結果。(a)根據(jù)人工目測統(tǒng)計結果作出的聚類圖;(b)本算法的聚類圖;(c)未進行Burst處理直接將所有鋒電位進行分類的結果;(d)根據(jù)ISI檢出Burst但不使用ICA分離處理而直接聚類的結果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)計了2 630個 Burst信號段，其中包含4類以上鋒電位的只占1.4%;因此，絕大多數(shù)候選Burst小段信號中只含有3類及以下鋒電位，滿足ICA源信號數(shù)量小于記錄信號通道數(shù)的限制條件。

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

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

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

圖4 候選Burst及其ICA分離示例。(a)信號中包含2類鋒電位和1個Burst;(b)信號中包含4類鋒電位和2個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ù)的聚類結果。(a)完全正確的理想分類圖;(b)本算法的分類圖;(c)未進行Burst處理直接將所有鋒電位進行分類的結果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ù)量和檢測分類結果(n=10)Tab.2 Results of the algorithm for 10 sets of synthetic signals

3 討論

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

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

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

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

(4)只對一部分信號小段進行ICA分析也大大減少了ICA計算量。長時間記錄信號分析時計算量過大也是阻礙ICA應用的一個大問題。因此，本研究的算法既發(fā)揮了ICA的優(yōu)勢又克服了其缺點。

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

4 結論

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

［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］ 封洲燕，光磊，鄭曉靜，等.應用線性硅電極陣列檢測海馬場電位和單細胞動作電位［J］.生物化學與生物物理進展，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ā)檢測自適應算法［J］.生物物理學報，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)元鋒電位檢測和分類的新方法［J］.生物化學與生物物理進展，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

國家自然科學基金(30770548，30970753).

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