基于模型分類算法的針刺神經電信號的統計分析

秦 晴，陳 偉，鄭海令

（天津職業技術師范大學天津市信息傳感與智能控制重點實驗室，天津300222）

針刺穴位能夠誘發神經系統產生響應活動，對神經和內分泌網絡進行調節，進而提高靶器官的各項功能[1-3]。在針刺動物實驗中，實驗者將電極植入動物脊髓，采集針刺神經響應信號[3-4]。將電極植入目標區，采集電極周圍神經元的胞外活動是神經科學中最常用的信號采集手段。由于電極采集的是周圍神經元集群的放電活動，所以實驗者的首要任務是對放電進行分類，將采集到的所有的動作電位分別歸入相應神經元[5]。盡管放電分類算法的研究困難且耗時，但仍是相當必要的，這是因為對于特定的輸入刺激，實驗者不能把采集到的所有電信號放在一起作為該刺激的響應，對于不同的輸入刺激，每個神經元個體不會呈現相同的放電活動。為了降低放電分類算法中的困難，早期實驗者采用比較精細的電極植入目標區，并放在給定神經元最近的位置，降低其他神經元的放電干擾，但是該方法存在一些缺陷：電極偏向于靠近高幅值放電，即放電明顯的神經元；由于電極的精細，其只能檢測一個或少數神經元的放電活動[5]。

隨著解剖學和生物物理學的快速發展，近期大多數的放電分類算法都是基于放電波形的。最早的基于放電波形的分類算法僅僅是基于動作電位振幅的分類處理，在很多情況下，振幅參數并不能將不同神經元的動作電位區分開。為此，基于放電波形的分類算法加入了時間窗的概念，定義1個或幾個時間窗為特定神經元的放電波形，然后在時間軸上進行滑動對比，形狀相同的放電波形就是該神經元的動作電位。但是，這種算法仍存在以下缺點：一是需要人為地選取作為參考的時間窗，這在成百個通道的分類處理中是不現實的；二是放電波形的疊加和背景噪聲都能夠導致放電波形發生變形，此時分類效果就會明顯變差。變異的放電波形被剔除，大量的放電信息被遺漏。為解決變異波形的分類問題，本文在波形分類算法的基礎上，采用一種基于模型的分類優化算法，該優化算法能夠有效地辨識變異的放電波形，降低動作電位的漏檢率，提高分類結果的準確性。

1 動作電位的分類算法

1.1 小波聚類算法

為了解決分類過程中人為導致的低效率性和誤差等問題[6]，文獻[7-8]提出了基于小波變換（wavelet transform）和超順磁聚類（super-paramagnetic clustering）的分類算法，簡稱小波聚類算法。算法分為3步：①設置閾值提取所有的動作電位；②對動作電位的放電波形進行小波變換，提取小波特征參數；③根據小波特征參數進行自動的超順磁聚類分析，獲取每個神經元的放電序列。

在第①步中，電極采集到的信號先經過一個高頻濾波器（300～3 000 Hz），然后根據濾波后的數據自動選取放電閾值，提取全部動作電位。閾值方程為

式中：x為經過高頻濾波的信號。

在第②步中，小波變換是對信號的時頻分析，以最佳的時間和最佳的頻域對信號做出最佳的解釋，同時它不需要信號時刻地保持穩定性。小波變換被定義為

其中，Ψa，b（t）是小波方程，小波變換是小波方程與信號x的卷積。根據方程（2）和（3）可以得到小波變換的逆變換，重構信號x。

在第③步中，超順磁聚類分析[7-8]基于概率統計知識的聚類方法，估計某個點與其臨近點的相互作用關系，這里的點yi表示動作電位i的m維小波特征。點間的互相作用強度被定義為

式中：‖yi-yj‖2為2個點的距離；K為臨近點的個數；a為臨近點的平均距離。

根據互相作用強度得到這2個點是同一類的概率

式中：T為溫度；〈δsi，s〉j為臨近點的狀態距離，決定了2個動作電位i和j是不是同一類，只有當〈δsi，s〉j≥0時，動作電位i和j才有可能是同一神經元的動作電位。

1.2 基于模型的優化算法

雖然小波聚類算法很好地解放了人為操作，極大提高了分類效率，但是在噪聲干擾[9-11]和波形疊加[12-14]問題的處理上并不是很理想。尤其在波形疊加的情況下，小波聚類方法的漏檢率很高。當神經元發生同步或近似同步放電時就會產生波形疊加，與背景噪聲相比，波形疊加并不是隨機事件，而是神經元集群網絡活動的重要體現。所以，由波形疊加所導致的分類錯誤是高度系統化的，不能被忽視。本文采用一種基于模型的分類優化算法，辨識同步或近似同步放電產生的疊加波形[15]。

假設t時刻電極記錄的動過電位是神經元集群中所有放電波形在t時刻以不同形式的線性疊加，建立t時刻動作電位的生成模型[15-16]為

式中：列向量v（t）=（vi（t））ne×1為t時刻所有電極記錄的動作電位；nc為多電極的個數；vi（t）為t時刻第i個電極記錄的動作電位；所有的放電波形都被平均地離散成nτ份，矩陣wj（τ）為被所有電極記錄的神經元j的第τ段波形；xj（t）為一個二進制變量，代表著神經元j在t時刻是否發生放電；η（t）為均值為0的背景白噪聲。

根據方程（6），得到整個觀測時間內的動作電位的生成模型

式中：W*X為方程（6）中的卷積運算；η（t）為一個多變量高斯白噪聲。

故V-W*X服從多變量高斯分布，那么在放電時刻和放電波形條件下，動作電位產生的條件概率函數為

式中：Λ是多維高斯變量η的協方差矩陣，它不僅表示時間上噪聲的協方差，還表示空間上電極采集位置所造成的噪聲的協方差。

將小波聚類算法的分類結果作為優化算法的初始值，定義每個神經元放電序列xj和放電波形wj（τ）的先驗概率函數。假設神經元放電是一個貝努利事件，那么xj（t）服從貝努利分布，且每個時間窗內都是獨立的，放電序列xj的先驗概率函數為

其中：xj（t）為一個二進制變量。t時刻神經元j產生放電時，xj（t）=1；否則xj（t）=0。從小波聚類方法的分類結果中可以獲取神經元j產生放電事件的先驗概率pj，如在針刺實驗中，動作電位采樣頻率為40 000 Hz。如果小波聚類算法的分類結果中，針刺神經元平均每秒產生10個動作電位，即平均4 000個采樣點就會產生1個放電，那么針刺神經元放電事件的先驗概率為1/4 000。

根據先驗函數、條件概率函數以及貝葉斯理論，動作電位由放電序列和放電波形誘發的后驗概率函數為

令該后驗概率函數取最大值，從而估計放電序列和放電波形[16－22]。

2 秩和檢驗

為了檢驗不同針刺手法的差異性，本文基于模型優化算法，對20次實驗的分類結果進行秩和檢驗。在秩和檢驗中，原假設H0是針刺補法和針刺瀉法所誘發的響應活動，沒有顯著差異。通過比較20次實驗中捻轉補法（“nb”）和捻轉瀉法（“nx”）所誘發的放電個數，可得到“nb”的秩和。同樣地，通過比較20次實驗中提插補法（“tb”）和提插瀉法（“tx”）所誘發的放電個數，可得到“tb”的秩和Ttb。由于樣本容量=20（nnb=nnx=ntb=ntx=20），因此秩和的分布接近正態分布。在這種情況下，可以采用Z檢驗，其基本公式為

若顯著性水平α定義為0.05，則根據正態分布的Z值表，可得H0得到的接受域為|Z|＜1.96。如果Znb和Ztb均不在該接受域內，H0被拒絕，這就意味著補法和瀉法所誘發的響應活動具有顯著差異性；否則H0是真，補法和瀉法無差異性。

3 針刺電信號的分類

3.1 2種分類方法的比較

采用針刺大鼠足三里實驗中脊髓背根處神經元集群的放電活動作為原始數據，分別進行2種算法的放電分類，并比較它們的分類效果。針刺試驗中，選擇健康的成年大鼠（體重180～200 g，年齡2～3歲）為實驗對象。實驗前，所有大鼠飼喂7 d，以適應實驗室標準環境。實驗過程中，先用20%氨基甲酸乙酯（1.5 g/kg）對大鼠進行深度麻醉，然后在解剖顯微鏡下，從L4背根的中樞端分離神經細束，將細束搭放在一對雙極鉑金絲記錄電極上，再用BIOPAC－MP150生理記錄儀對針刺誘發的脊髓背根神經元的放電活動進行記錄。確定足三里穴所對應的脊髓背根神經纖維（L4）后，待神經纖維靜息后，依次對足三里穴進行不同手法的針刺刺激。針刺手法分別為捻轉補法（“nb”），捻轉瀉法（“nx”），提插補法（“tb”）和提插瀉法（“tx”），針刺頻率為100次/min。為排除針刺順序對實驗結果產生的影響，在實驗過程中采取隨機的施針順序。每種手法刺激時間為0.5 min，然后留針5 min，待大鼠靜息后，再進行另一種手法的針刺刺激。

小波聚類算法的結果如圖1所示。圖1（a）—（d）分別為“nb”、“nx”、“tb”和“tx”4種針刺手法誘發的原始數據。4種手法條件下所有放電波形的小波特征參數，如圖1（e）所示，這是一個二維的散點圖。圖中小波特征參數分布在4個不同顏色的區域，所以這些動作電位相應地被分為4類，每一類的平均動作電位如圖1（f）所示。在此，定義紅色高幅值的動作電位是由神經元1產生，綠色余弦動作電位是由神經元2產生，藍色正弦動作電位是由神經元3產生，青色低幅值放電是由神經元4產生。由于神經元1放電波形的幅值明顯較高，所以神經元1的小波特征參數的分布區域離其他3類比較遠。4種手法條件下，估計得到的4個神經元的放電率如表1所示。然后，將小波聚類算法的分類結果作為初始值，根據模型算法，對4組原始數據進行重新分類，估計結果如表2所示。

表1 小波聚類算法中4個神經元在4種針刺手法下的放電率 個/s

表2 模型算法中4個神經元在4種針刺手法下的放電率 個/s

圖1 小波聚類算法的結果

從表1中可以看出，在不同針刺手法條件下4個神經元放電率的差異性不大。從表2可知，在模型算法中除了神經元4，其他3個神經元被檢測到的放電率均顯著增加，有的甚至是小波聚類算法中放電率的2倍。而對于第4類動作電位，由于其幅值最低，且基電位過高，所以被認為不是神經元的響應活動，而是周圍噪聲的結果，所以4類動作電位不予考慮。為了更直觀地體現模型算法的優點，選取2種算法在一次針灸刺激（約0.6 min）中的分類結果進行比較，電極采集的原始數據如圖2所示，模型算法檢測到的4個神經元的放電序列如圖3所示，小波聚類算法檢測到的4個神經元的放電序列如圖4所示。

圖2 電極采集的原始數據

從圖3和圖4可知，藍色圓點是2種算法都檢測到的動作電位。在圖3中，紅色圓點是模型算法多檢測的動作電位，模型算法很好地解決了多個神經元放電疊加的問題。尤其是紅色虛線框的部分，原始數據的膜電壓集中地大幅度波動，表明神經元集群的響應活動強烈，此時模型算法能將所有的動作電位都檢測出來，而基于放電波形的小波聚類算法卻因為放電波形的多次疊加沒有辨識出一個動作電位。另外，盡管小波聚類算法已經很好地解決了噪聲的問題，但是與模型算法相比，仍有不足，如8.7 s附近的紅色圓點（神經元1的高幅值放電），8.9～9 s和9.2 s附近的紅色圓點（神經元3的放電），小波聚類算法都沒有檢測到，其所對應的膜電壓曲線都存在動作電位的明顯波動。在圖4中，綠色圓點是模型算法漏掉的第4類動作電位。對于模型算法漏掉的動作電位（綠色圓點），其絕大部分分布在針灸刺激結束后的時間段，如綠色虛線框，同時其所對應的膜電壓曲線在背景噪聲范圍內，所以這類動作電位被剔除。

圖3 模型算法檢測到的4個神經元的放電序列

圖4 小波聚類算法檢測到的4個神經元的放電序列

3.2 4種針灸手法條件下神經元放電事件的統計分析

為了更加深入地研究不同手法刺激條件下神經元響應的差異性，本文運用模型分類算法對20次實驗的原始數據進行分類處理，并對分類結果進行秩和檢驗，比較針刺補法和針刺瀉法的差異性。秩和檢驗結果如表3和表4所示。

表3 “nb”和“nx”差異性的秩和檢驗

表4 “tb”和“tx”差異性的秩和檢驗

從表3和表4可以看出，捻轉補法和提插補法條件下，3個神經元的Z值不在接受域內，均拒絕。這就意味著捻轉補法和捻轉瀉法條件下3個神經元的放電活動都具有明顯的差異性，提插補法和提插瀉法條件下3個神經元的放電活動也具有明顯差異性。

為了進一步揭示補法和瀉法差異性的具體體現，本文對20次實驗中神經元集群的放電事件進行統計分析。4種針刺手法條件下，3個神經元在20次實驗中放電事件的統計結果如表5所示。

表5 4種針刺手法條件下3個神經元在20次實驗中放電事件的統計結果 個/s

從表5可知，補法條件下3個神經元20次實驗的平均放電率都分別大于相應瀉法條件下3個神經元放電的平均放電率，其中第5行表示捻轉補法條件下3個神經元在單位時間內多產生的放電個數，第6行表示提插補法條件下單位時間內多產生的放電個數。由此可得，補法和瀉法條件下神經元響應的差異性在于補法能夠誘發神經元產生較多的放電響應活動。

4 結語

本文根據針刺實驗中電極記錄的原始信號，采用基于模型的分類優化算法，克服了傳統分類算法中波形變異和波形疊加問題，將脊髓背根處3個神經元的放電波形和放電序列準確高效地估計出來。同時，對多次實驗中不同針刺手法條件下神經元放電事件進行秩和檢驗和統計分析。結果發現，針刺補法和針刺瀉法產生的神經響應活動存在明顯的差異性，在針刺補法條件下神經元放電響應活動明顯多于針刺瀉法。