基于局部場電位固有模態(tài)分量的響應(yīng)調(diào)諧特性研究

萬 紅 張曉娜 劉新玉 李曉燕

(鄭州大學(xué)電氣工程學(xué)院，鄭州 450001)

引言

調(diào)諧特性是初級視覺皮層(primary visual cortex，又稱V1 區(qū))神經(jīng)元感受野的主要特性之一［1］。傳統(tǒng)的對V1 區(qū)神經(jīng)元對空間頻率響應(yīng)特性的研究，大都基于神經(jīng)元鋒電位(spike)發(fā)放率來獲得神經(jīng)元對刺激響應(yīng)的調(diào)諧曲線。然而在微電極陣列胞外采集技術(shù)中，spike 受到噪聲的干擾經(jīng)常會出現(xiàn)漏檢和誤檢，即假陰性和假陽性問題，使后續(xù)分析受到嚴重影響。因此，與spike 活動同時記錄到的局部場電位信號，逐漸成為目前視覺信息處理機制研究的熱點之一［2］。

局部場電位(local field potential，LFP)是記錄電極尖端附近局部區(qū)域的神經(jīng)元興奮性和抑制性突觸后電位的總和。近年來，LFP 對于視覺刺激特征的調(diào)諧作用受到廣泛關(guān)注。Gawne 等指出，局部場電位信號與鋒電位發(fā)放率對于視覺刺激有著相似的調(diào)制作用，當(dāng)鋒電位發(fā)放率達到最大時，局部場電位的能量同時也達到了峰值［3］;Ince 等通過對朝向、空間頻率、時間頻率等3 個視覺刺激特征的研究，表明V1 區(qū)LFP 頻帶能量能夠?qū)ι鲜鎏卣鬟M行有效編碼［4］;Belitski 等也發(fā)現(xiàn)LFP 信號不同頻帶編碼了不同的視覺刺激特征［5］。諸多研究表明LFP對不同視覺刺激特征具有調(diào)諧功能，然而對于LFP中與刺激的最佳響應(yīng)范圍，研究者們沒有達成共識。

Henrie 等通過對麻醉猴的研究發(fā)現(xiàn)，LFP 在刺激前后功率譜會發(fā)生明顯的變化，并且不同頻帶的能量變化不同，能量變化最明顯頻帶主要集中在25～90 Hz［6］;Kavser 等以清醒貓為對象，研究了LFP不同頻段對時間頻率、空間頻率和朝向的選擇性，結(jié)果表明具有調(diào)諧特征的頻帶，主要集中在8 ～23 Hz 和39 ～109 Hz［7］;Philipp 等通過記錄清醒獼猴主視區(qū)的LFP 信號，研究了神經(jīng)元感受野的特征選擇性，發(fā)現(xiàn)LFP 最佳調(diào)諧范圍主要集中在30 ～90 Hz［8］。Magri 等通過記錄麻醉獼猴初級視覺皮層的信號，將LFP 劃分為表征自然視覺刺激的多個頻帶［9］。就目前研究結(jié)果來看，由于動物類型差異和LFP 的非平穩(wěn)性，使得不同的實驗獲得的最佳響應(yīng)范圍并不一致。因此，獲得LFP 最佳響應(yīng)頻帶范圍、提取出神經(jīng)元響應(yīng)特征，還需要深入研究。

局部場電位是典型的非線性、非平穩(wěn)信號，經(jīng)典的時頻方法處理效果并不理想。近些年來，具有自適應(yīng)特性的希爾伯特黃變換(Hilbert-Huang transform，HHT)正逐漸得到研究者的青睞。HHT 沒有先驗的基函數(shù)，而且分解過程完全基于數(shù)據(jù)本身［10］。Liang 等的研究表明，HHT 是神經(jīng)電信號分析的強有力工具［11-13］。本研究基于HHT 算法，首先將LFP 信號進行經(jīng)驗?zāi)B(tài)分解，將其第二階固有模態(tài)分量進行Hilbert 變換并求取其能量，利用該能量特征研究LFP 對不同空間頻率刺激的調(diào)諧特性，并與同一電極記錄到的多神經(jīng)元鋒電位(multi-unit activity，MUA)和小波分解提取的Gamma 頻帶進行了對比。

1 材料和方法

1.1 材料

實驗動物選用年齡12 周左右、體重200 ～250 g的6 只成年Long Evans(LE)大鼠，雌雄不拘。術(shù)前腹腔注射10%水合氯醛麻醉，實施開顱手術(shù)。參考大鼠腦立體定位圖譜，確定初級視覺皮層的中心位置，用顱鉆進行開顱手術(shù)。用顯微操縱器(MX7600，Siskiyou Inc.)植入2 ×8 微電極陣列(Microprobe，鉑銥合金，電極間距250 μm，電極直徑125 μm，尖端直徑15 μm，電極阻抗在0.5 ～1.0 MΩ 之間)，并用數(shù)據(jù)采集系統(tǒng)(CerebusTM，Blockrock Inc.)記錄大鼠V1 區(qū)神經(jīng)元誘發(fā)放電信號，包括檢測到的LFP 和MUA 信號。MUA 將單個通道的spike 作為一個整體，是電極尖端附近多個神經(jīng)元胞體的活動［14］。視覺刺激采用不同空間頻率的全屏運動光柵，空間頻率分別取0.01、0.02、0.05、0.1、0.2、0.4 Cpd，時間頻率為2 Hz，刺激1 s，間隔1 s。每組刺激各空間頻率光柵隨機出現(xiàn)，每次實驗重復(fù)20 組。

1.2 方法

1.2.1 希爾伯特黃變換

HHT 包括經(jīng)驗?zāi)B(tài)分解(empirical mode decomposition，EMD)和Hilbert 變換兩個部分。首先采用EMD 方法將LFP 信號分解為若干固有模態(tài)分量(intrinsic mode function，IMF)，然后對每個IMF 進行Hilbert 變換，求出其瞬時頻率和幅值，構(gòu)成被分析信號的時間-頻率-幅值(能量)的三維譜圖，即Hilbert 譜。EMD 分解的IMF 須滿足以下兩個條件:一是在整個信號長度上，極值點和過零點數(shù)必須相等或至多相差一個;二是在任何一點，由局部極大值點形成的包絡(luò)線和局部極小值點形成的包絡(luò)線的均值為零。

LFP 的HHT 變換過程如下:

首先，對其中一個通道LFP 信號x 進行EMD 分解，得到從高頻到低頻的n 個IMF 分量c1，c2，…，cn和余項rn，即

然后對每個IMF 分量進行Hilbert 變換

當(dāng)y(t)與c(t)形成一個復(fù)共軛時，就可構(gòu)造解析信號z(t)

式中，α 代表瞬時幅值，θ 代表瞬時相位，并且有

從而得到IMF 分量的瞬時頻率ωi(t)

式(4)和式(5)表明信號的幅值和瞬時頻率都是時間的函數(shù)，則原始信號可表示為

取其實部，即原始信號的Hilbert 譜，

式中，Re 表示實部，這里省略了殘余函數(shù)rn，因為它是一個單調(diào)函數(shù)或者常量，代表著長期振動。

1.2.2 響應(yīng)特征提取

瞬時能量(instantaneous energy，IE)提供了一種能量隨時間變化的信息，由Hilbert 譜可以得到

式中，H 為Hilbert 變換，ω 表示目標(biāo)頻帶，t 代表時間。通過式(8)可以計算單一模態(tài)的瞬時能量。

1.2.3 調(diào)諧指數(shù)

調(diào)諧指數(shù)度量了信號特征對刺激參數(shù)的響應(yīng)程度。為了定量的表示信號的調(diào)諧性能，分別計算了MUA 和LFP 在不同空間頻率刺激下的調(diào)諧指數(shù)。MUA 的調(diào)諧指數(shù)dMUA為［8］

式中，Pmax表示MUA 發(fā)放率調(diào)諧曲線的最大值，Pmin表示MUA 發(fā)放率調(diào)諧曲線的最小值。同理，LFP 的調(diào)諧指數(shù)dLFP為［8］

式中，Emax表示LFP 能量調(diào)諧曲線的最大值，Emin表示LFP 能量調(diào)諧曲線的最小值。

為了定量衡量不同方法獲得的神經(jīng)元響應(yīng)特征的一致程度，引入了神經(jīng)元響應(yīng)一致率指標(biāo)，定義為:設(shè)總神經(jīng)元個數(shù)為M，其中兩種信號對最佳刺激空間頻率響應(yīng)一致的神經(jīng)元個數(shù)為N，那么最佳空間頻率刺激響應(yīng)的一致率λ 定義為

1.2.4 小波分解提取LFP 的Gamma 頻帶

小波分解結(jié)果依賴于所選擇的小波基。基于LFP 的特點，選擇具有緊支正交基，滿足精確信號重建條件的“db5”作為小波基函數(shù)［15］。小波分解表示為

小波重構(gòu)

式中，小波基函數(shù)為

式中，a0＞0，j，k ∈Z，f(t)為原始信號，為尺度因子，kτ0為平移因子。

為了提取LFP 的Gamma 頻帶，通過式(13)將LFP 信號進行5 層小波分解，將其中的62.5 ～125 Hz 頻帶再進行分解得到62.5 ～93.75 Hz，與第一次分解得到的31.25 ～62.5 Hz 進行小波重構(gòu)，從而得到Gamma 頻帶(30 ～90 Hz)信號。

2 結(jié)果

2.1 各階固有模態(tài)分量特征分析

首先對記錄到的16 通道原始數(shù)據(jù)進行低通濾波，通帶截止頻率200 Hz，階數(shù)200，利用零均值處理去除基線漂移，記為LFP。預(yù)處理前后的LFP 信號如圖1 中(a)和(b)所示，其中刺激時間為1 s。將除余量之外的各階IMF 信號合成求取Hilbert 譜，如圖1(c)所示。由圖可知，LFP 在刺激前后能量發(fā)生明顯變化(t-test，P ＜0.05)。

為了進一步提取LFP 的響應(yīng)頻帶范圍，對LFP進行了EMD 分解，得到7 個固有模態(tài)分量和一個余量，如圖1(d)所示。從分解結(jié)果來看，每一階分量都具有不同振幅和頻率，第一階模態(tài)分量的頻率最高，隨著模態(tài)階數(shù)的增加其所含的周期變長，直到最后一階逼近一條線性的時間函數(shù)，而且前幾個分量包含了信號的主要特征，在刺激后能量有不同程度的增加。

2.2 IMF2 的調(diào)諧性能

調(diào)諧特性的傳統(tǒng)研究方法，大都基于spike 發(fā)放率，而且已得到了大量文獻的證實［16-17］。因此選用MUA 發(fā)放率為基準(zhǔn)，來統(tǒng)計分析各階固有模態(tài)分量是否具有調(diào)諧特性。表1 給出了基于各階固有模態(tài)分量，對不同空間頻率的調(diào)諧指數(shù)平均值及與MUA 的一致率。從表中可以發(fā)現(xiàn)，LFP 的第二階固有模態(tài)分量(IMF2)的調(diào)諧指數(shù)最大，其與MUA 最佳空間頻率調(diào)諧的一致率為68.75%，顯著高于其它固有模態(tài)分量。這一結(jié)果表明，LFP 攜帶的刺激相關(guān)信息主要集中在IMF2。

以IMF2 作為特征頻帶，以式(8)定義的計算單一模態(tài)的瞬時能量為能量特征，分析了大鼠V1 區(qū)空間頻率響應(yīng)特性，并將其與MUA 發(fā)放率進行了對比，結(jié)果如圖2 所示。

圖1 光柵刺激前后LFP 響應(yīng)特性。(a)原始LFP 信號(刺激時間1 s);(b)預(yù)處理后LFP 信號;(c)LFP 的Hilbert譜;(d)LFP 經(jīng)驗?zāi)B(tài)分解后的各階固有模態(tài)分量(C1 －C7)和余量R8Fig．1 The response characteristics of LFP before and after grating stimulation． (a)The original LFP signal(stimulation time is one second);(b)LFP signal after pretreatment;(c)The Hilbert spectrum of LFP;(d)The intrinsic mode functions (C1 －C7)and residual vector R8 of LFP after empirical mode decomposition

表1 LFP 固有模態(tài)分量調(diào)諧指數(shù)統(tǒng)計表Tab．1 The statistical table of tuning index based on the intrinsic mode functions of LFP

圖2(a)表示20 次重復(fù)刺激IMF2 和MUA 發(fā)放率的對比結(jié)果。由圖可知，IMF2 幅值和MUA 的發(fā)放率在刺激前后都發(fā)生了顯著變化，兩者具有一致性。圖2(b)是20 次重復(fù)刺激IMF2 在灰屏和光柵刺激能量的對比圖，虛線表示灰屏?xí)r的能量變化，實線表示光柵刺激時的能量變化，陰影部分為標(biāo)準(zhǔn)差，顯著水平為0.05。由圖可知光柵刺激下IMF2能量明顯增加，進一步表明IMF2 攜帶了與刺激相關(guān)的信息。圖2(c)為IMF2 和MUA 在不同空間頻率刺激下的標(biāo)準(zhǔn)化調(diào)諧曲線(z-score 標(biāo)準(zhǔn)化，Matlab)，二者對不同空間頻率刺激都有不同的響應(yīng)強度，而且調(diào)諧曲線具有較好的一致性，同時在空間頻率為0.02 Cpd 時響應(yīng)最強烈。

2.3 IMF2 功率譜估計

為確定IMF2 的頻帶范圍，對其進行了功率譜估計，并與小波分解提取的30 ～90 Hz 的Gamma 頻帶的功率譜進行了對比，如圖3 所示。圖中選用Welch 方法分別對IMF2 和Gamma 進行功率譜估計，Welch 估計的關(guān)鍵是選取合適的窗函數(shù)和窗長，選取主瓣較窄的Hamming 窗［15］。經(jīng)對比研究發(fā)現(xiàn):選用Hamming 窗分辨率20 Hz，即對應(yīng)窗長N =200，能夠準(zhǔn)確提取頻帶能量。由圖可知IMF2 的頻帶范圍主要集中在35 ～80 Hz 內(nèi)，比Gamma 頻帶范圍更窄。

2.4 IMF2 與MUA、小波分解的調(diào)諧性能對比

圖2 IMF2 與MUA 空間頻率調(diào)諧特性對比。(a)20 次重復(fù)刺激下IMF2 和MUA 發(fā)放率;(b)灰屏和光柵刺激下IMF2 能量變化;(c)IMF2 與MUA 調(diào)諧曲線Fig．2 Comparison of spatial frequency tuning characteristics between IMF2 and MUA． (a)IMF2 and MUA firing rate under 20 times repeat stimulation;(b)The change of IMF2 energy under gray screen and grating stimuli;(c)Tuning curves of IMF2 and MUA

圖3 IMF2 和Gamma 頻帶的功率譜估計Fig．3 The power spectrum of IMF2 and Gamma-band

為進一步驗證該方法在提取特征頻帶時的優(yōu)越性，利用從6 只LE 大鼠記錄到的共63 通道神經(jīng)元響應(yīng)信號，計算了IMF2 與MUA、小波分解提取的Gamma 頻帶對空間頻率刺激響應(yīng)具有一致性的通道調(diào)諧指數(shù)，結(jié)果如圖4 所示。圖4(a)為IMF2 和MUA 的調(diào)諧指數(shù)對比圖。圖中橫坐標(biāo)和縱坐標(biāo)分別為IMF2 和MUA 的調(diào)諧指數(shù)，十字表示平均值。由圖可知，IMF2 的平均調(diào)諧指數(shù)為0.795 1，MUA為0.631 3，這一結(jié)果表明IMF2 調(diào)諧性能要強于MUA。圖4(b)是IMF2 與小波分解提取的Gamma頻帶的調(diào)諧指數(shù)對比圖。圖中橫坐標(biāo)和縱坐標(biāo)分別為IMF2 和小波分解提取的Gamma 頻帶的調(diào)諧指數(shù)，十字表示平均值。由圖可知，前者平均調(diào)諧指數(shù)為0.795 1，后者為0.664 6。由此可見，基于HHT方法提取的LFP 響應(yīng)頻帶較小波分解方法更為有效。

3 討論

LFP 活動反映了大腦皮層局部區(qū)域神經(jīng)元集群的響應(yīng)特性，可以從更精細的時間尺度和空間尺度上解釋認知的過程，在視覺信息處理機制研究中得到了廣泛應(yīng)用。近年來研究發(fā)現(xiàn)，LFP 中與刺激相關(guān)的響應(yīng)頻帶主要集中在Gamma 帶，然而由于LFP的非平穩(wěn)性，使得不同文獻報道的響應(yīng)頻帶范圍并不一致。因此，確定LFP 的響應(yīng)范圍、提取LFP 的響應(yīng)特征，還需要進一步深入研究。

圖4 不同方法得到的調(diào)諧指數(shù)對比圖。(a)IMF2 與MUA 調(diào)諧指數(shù)對比圖;(b)IMF2 與小波分解的Gamma 頻帶調(diào)諧指數(shù)對比圖Fig．4 Comparison of the tuning index obtained from different ways．(a)Comparison of the tuning index of IMF2 and MUA;(b)Comparison of the tuning index of IMF2 and Gamma-band extracted with wavelet decomposition

本研究基于HHT 變換，研究了LFP 的響應(yīng)頻帶范圍。由于HHT 變換引入了EMD 分解，使得經(jīng)EMD 分解得到的各階固有模態(tài)分量完全基于信號本身的特性，無需先驗基函數(shù)，可以有效反應(yīng)信號的固有特征，其分解的客觀性、自適應(yīng)性凸顯了該方法在非平穩(wěn)信號處理中的優(yōu)勢，具有廣闊的應(yīng)用前景。由實驗結(jié)果可以看出，HHT 能夠有效分解出LFP 的響應(yīng)頻帶(如圖1)，彌補了傳統(tǒng)小波算法需要選取合適小波基的不足，拓展了LFP 處理的新思路，為LFP 信號的響應(yīng)特征提取及后續(xù)研究奠定了基礎(chǔ)。

本研究得出的結(jié)論與現(xiàn)有的研究結(jié)果比較吻合。目前的研究表明，LFP 的Gamma 頻帶(30 ～90 Hz)能量能夠表征視覺刺激特征，本研究表明LFP經(jīng)過EMD 分解得到的IMF2 的頻帶范圍主要集中在35 ～70 Hz，這個范圍的主要頻帶與Gamma 帶重疊，而且較Gamma 頻帶范圍更窄、更精確(如圖3)，并且其調(diào)諧指數(shù)高于Gamma 頻帶(如圖4(b))。這些結(jié)果表明基于HHT 方法提取LFP 響應(yīng)信號更為有效。激的調(diào)諧特性，并與同一電極記錄到的MUA 和小波分解提取的Gamma 頻帶進行了對比。

通過對記錄的63 個通道神經(jīng)元統(tǒng)計分析，LFP的IMF2 攜帶了刺激信息，與其他模態(tài)分量相比，IMF2 對刺激光柵空間頻率的調(diào)諧特性最強，與MUA 空間頻率響應(yīng)的一致率達68.75%，且其調(diào)諧指數(shù)的均值(0.795 1)高于 MUA (0.631 3)和Gamma 頻帶(0.664 6)。表明基于LFP 第二固有模態(tài)分量來研究神經(jīng)元感受野的最佳空間頻率調(diào)諧特性，比MUA 和小波提取的Gamma 頻帶更有優(yōu)勢，為進一步從LFP 分析視覺皮層信息處理機制，提供了另一種可行的有效方法。

4 結(jié)論

本研究基于HHT 算法，首先將LFP 信號進行經(jīng)驗?zāi)B(tài)分解，然后將其第二階固有模態(tài)分量作為特征頻段，進行Hilbert 變換并求取其能量特征，研究大鼠初級視覺皮層LFP 對不同空間頻率光柵刺

［1］ Mazer JA，Vinje WE，McDermott J，et al． Spatial frequency and orientation tuning dynamics in area V1 ［J］. Neurobiology，2002，99(3):1645 -1650.

［2］ Belitski A，Panzeri S，Magri C，et al． Sensory information in local field potentials and spikes from visual and auditory cortices:time scales and frequency bands ［J］. Journal Computational Neuroscience，2010，29(3):533 -545.

［3］ Gawne TJ. The local and nonlocal components of the local field potential in awake primate visual cortex ［J］. Journal of Computational Neuroscience，2010，29(3):615 -623.

［4］ Ince RAA，Mazzoni A，Bartels A，et al. A novel test to determine the significance of neural selectivity to single and multiple potentially correlated stimulus features［J］. Journal of Neuroscience Methods，2012，210(1):49 -65.

［5］ Belitski A，Gretton A，Magri C. Low-frequency local field potentials and spikes in primary visual cortex convey independent visual information［J］. The Journal of Neuroscience，2008，28(22):5696 -5709.

［6］ Henrie JA，Shapley R. LFP power spectra in V1 cortex:the graded effect of stimulus contrast ［J］. Journal of Neurophysiology，2005，94(1):479 -490.

［7］ Kayser C，Konig P. Stimulus locking and feature selectivity prevail in complementary frequency ranges of V1 local field potentials［J］. Neuroscience，2004，19(2):485 -489.

［8］ Philipp B，Georgios A. Comparing the feature selectivity of the gamma-band of the local field potential and the underlying spiking activity in primate visual cortex ［J］. Frontiers in Systems Neuroscience，2008，2(2):1 -11.

［9］ Magri C，Mazzoni A，Logothetis NK，et al． Optimal band separation of extracellular field potentials ［J］. Tournal of Neuroscience Methods，2012，210(1)，66 -78.

［10］ Huang NE，Shen Z，Long SR，et al． The empirical mode decomposition and the hilbert spectrum for nonlinear and nonstation time series analysis［J］. Proceedings of the Royal Society of London，1998，454:903 -995.

［11］ Liang HL，Steven LB，Desimone R，et al． Empirical mode decomposition:a method foranalyzing neural data ［J］.Neurocomputing，2005，65 －66:801 -807.

［12］ Sweeney CM，Nasuto SJ. A novel approach to the detection of synchronisation in EEG based on empirical mode decomposition［J］.Journal of Computational Neuroscience，2007，23(1):79 -111.

［13］ 袁玲，楊幫華，馬世偉. 基于HHT 和SVM 的運動想象腦電識別［J］. 儀器儀表學(xué)報，2010，31(3):649 -654.

［14］ Kajikawa Y，Schroeder CE. How local is the local field potential?［J］. Neuron，2011，72(5):847 -858.

［15］ 尚志剛，馮平艷，劉新玉等. 局部場電位γ 頻帶能量對朝向調(diào)諧特性研究［J］. 鄭州大學(xué)學(xué)報(工學(xué)版)，2012，33(6):28 -32.

［16］ Molotchnikoff S，Gillet PC. Spatial frequency characteristics of nearby neurons in cat's visual cortex［J］. Neuroscience Letters.2007，418(3):242 -247.

［17］ Cristopher M，Tryker MP. Highly selective receptive fields in mouse visual cortex［J］. The Journal of Neuroscience，2008，30(28):7520 -7536.