獼猴伸展抓握過程中的spike與LFP鎖相相位分布特征分析

劉秀玲，郭天翔，王潔，孫洪吉，熊鵬

(1. 河北大學 電子信息工程學院，河北 保定 071002；2. 河北省數字醫療工程重點實驗室，河北 保定 071002；3.軍事科學院軍事醫學研究院 軍事認知與腦科學研究所，北京 100850)

大腦神經元的主要活動包括單神經元的動作電位和由突觸電位形成的局部場電位.單神經元動作電位反映單細胞的電活動，這種通過神經軸突流動的電信號是神經信息傳輸的基礎；局部場電位反映一定范圍內的局部電活動.相位是LFP(localfield potential)的重要特征之一，鎖相即相位鎖定，表示2個信號間的相位同步性.spike的發放和LFP的節律之間的鎖相關系包含了重要的神經編碼信息[1].

Alex等[2]利用spike與LFP的beta節律的鎖相特點對獼猴行為狀態和腦網絡活動變化性的關系進行了分析. García等[3]利用spike和LFP的鎖相研究在聽覺皮層對刺激的處理中spike和LFP的相干性動力學. James等[4]利用2個腦區特定功能相關頻帶下多神經簇的發放與LFP的鎖相關系以研究大鼠癲癇模型中梨狀皮層和內側丘腦的放電模式. Sritharan等[5]將spike和場電位的相干性和鎖相統計數據作為量化研究初級感覺皮層神經元對2種不同紡錘波的作用. Guillem等[6]在研究大鼠前額葉皮層可卡因的使用與非藥物獎勵的個體偏好中，使用前額葉皮層的單神經元活動與theta段振蕩的鎖相時間機構作為研究方法之一.馬維建等[7]在大鼠海馬體中用鎖相關系的變化證實了深部腦刺激技術所使用的電脈沖高頻刺激可以改變spike的病理性發放.

因此，spike與LFP的鎖相關系的研究對理解大腦神經信息編碼具有重要意義.以上對鎖相關系的研究證實了在多種大腦活動中spike與LFP鎖相的存在性，并且利用鎖相關系對大腦活動做了進一步分析，但關于鎖相關系中鎖相相位的分布特征的分析還比較匱乏.

基于這種情況，本文利用獼猴伸展抓握運動范式，在多個腦區內，對獼猴在實際運動中的峰電位與局部場電位的鎖相相位變化進行了分析，為利用鎖相解析運動中的神經編碼信息提供了新思路.

1 材料

1.1 實驗信號采集

本實驗數據來自2只獼猴：獼猴A和獼猴B，M1、S1、PPC共3個腦區的3 d實驗的spike發放時間數據(spike timestamp)和LFP信號，所有實驗動物和實驗數據均來源于軍事科學院軍事醫學研究院實驗動物中心[8].

剔除噪聲過大數據后，共使用獼猴A 77個通道的數據(M1 31通道,S1 21通道，PPC 25通道)，使用獼猴B 78個通道的數據(M1 31通道，S1 31通道，PPC 16通道)[8].

1.2 實驗范式

獼猴的實驗范式按照圖1所示流程進行.

圖1 實驗范式流程Fig.1 Process of experimental paradigm

如圖2，獼猴按照提示燈1和提示燈2的指示進行實驗，實驗范式分為運動準備階段(wait)、反應階段(reaction)、伸手階段(reach)和抓握階段(grasp).各階段之間由事件CenterHit、TargetOn、CenterRelease、TargetHit進行劃分.主要流程為：

1)水平面板上，提示燈1點亮，提示獼猴可以將右手放置在觸摸板上.當獼猴放置后，觸摸板感受到獼猴手部壓力，則為觸摸成功，此時的事件命名為CenterHit事件，運動進入等待階段；

2)當獼猴將手放在觸摸板上維持500 ms后，水平面板上的提示燈1熄滅，位于垂直面板上的提示燈2點亮，提示獼猴可以開始抓握物體，事件命名為TargetOn事件，此時等待階段結束，運動進入反應階段；

圖2 獼猴實驗平臺Fig.2 Rhesus monkey experiment platform

3)在獼猴接收到提示燈2點亮的提示后，手部離開觸摸板需要一定的反應時間，當獼猴手部離開觸摸板，伸向抓握物體時，觸摸板感受不到壓力，此時事件為CenterRelease事件，此時反應階段結束，運動進入伸手階段；

4)在獼猴手部觸摸到目標物體時，此時事件為TargetHit事件，此時伸手階段結束，運動進入抓握階段；

5)獼猴手部觸摸到目標物體，并保持抓握姿勢200 ms，則認為抓握運動完成，此時抓握階段結束，一次完整實驗成功[8].

2 方法

2.1 鎖相相位的計算

本文使用MATLAB對LFP信號進行了delta段(1～4 Hz)的濾波，并按照事件CenterHit和TargetHit+0.2 s的時間標簽(timestamp)，對全部的spike timestamp和LFP信號進行運動時間范圍內的截取，即保留

TCenterHit

(1)

其中X表示spike timestamp或LFP的數據，TCenterHit和TTargetHit表示事件CenterHit和TargetHit發生時的timestamp.

完成后對LFP信號進行Hilbert變換[9]提取瞬時相位，具體計算過程如下：對于LFP信號y(t),y表示時間t點上的場電位值，其定義解析信號z(t)表示為

z(t)=y(t)+iHy(t),

(2)

其中虛部Hy(t)是y(t)的希爾伯特變換，則信號的瞬時相位為

φ(t)=Atg(z(t)).

(3)

為方便計算，本文使用CircStat工具箱將所有的角度轉化為弧度[10]，角度轉化為弧度的計算公式為

(4)

其中αdeg表示相位角度值，αrad表示相位弧度值.對小于0的弧度值加上2π后調整為正數，保證所有的LFP相位弧度值在0～2π內，便于后續分析.接下來按照通道為所有spike timestamp找到對應的最近時間點的LFP delta節律的相位，若找到的時間點與spike timestamp的時刻相差超過0.01 s，則認為該spike timestamp無對應的鎖相相位.將所有的spike timestamp對應的鎖相相位按照腦區和通道分類進行保存用于后續的分析.

2.2 不同運動階段鎖相相位分布

按照CenterHit、TargetOn、CenterRelease、TargetHit 4個事件的timestamp對鎖相相位進行運動階段的劃分.

TCenterHitTTargetOnTCenterReleaseTTargetHit

(5)

其中，Xwait、Xreaction、Xreach和Xgrasp表示截取后等待、反應、伸手、抓握4個運動階段的數據，TCenterHit、TTargetOn、TCenterRelease、TTargetHit表示CenterHit、TargetOn、CenterRelease、TargetHit 4個事件的timestamp.

每個腦區內所有通道的鎖相相位疊加起來即代表該腦區的所有鎖相相位，使用MATLAB軟件函數ksdensity對各運動階段的鎖相相位進行概率密度分布計算以分析其分布特征，分析范圍為0～2π，指定評估點為0～2π內的步長為0.02π的所有相位點；各運動階段的相位通過Wilcoxon秩和檢驗進行兩兩階段之間的統計學檢驗判斷其是否具有顯著差異性，當P<0.05時，認為其具有顯著性差異[11].

2.3 事件前后100 ms內鎖相相位分布

按照TargetOn、CenterRelease、TargetHit 3個事件的timestamp±100 ms內截取3個腦區的鎖相相位，并對比分析2個窗口內的鎖相相位概率密度分布的特點和峰值位置變化.對2個窗口內的鎖相相位，基于其數據量較小的特點進行置換檢驗分析，判斷兩者之間是否具有顯著性差異，當P<0.05時，認為具有顯著性差異[12].

2.4 連續時間窗對鎖相相位概率密度分布的分析

為了進一步分析鎖相相位隨時間的變化情況，使用窗口長度(binlength)為100 ms，步長為25 ms的連續時間窗對獼猴連續運動過程中的鎖相相位變化進行分析，窗口的起點(binstart)與CenterHit對齊，為保證數據整齊，窗口終點(binend)應該為binstart加上平均實驗時長(triallength)，即平均實驗時長為所有單次范式完成時長的平均值，則每個窗口內的數據為

bini

(6)

其中，i表示窗口的序數，xi表示當前窗口內的數據，bini表示當前窗口的起始點，當i為1時，bini=binstart，同時bini的值也代表當前窗口的時間點，當bini

對每個時間窗口內的鎖相相位取均值作為該時刻鎖相相位，即Xi=mean(xi)，Xi表示當前窗口對應時刻的相位平均值.將所有通道的時間窗結果疊加作為所在腦區的時間窗分析結果，對每個窗口內的鎖相相位進行概率密度分布分析，分析其鎖相相位分布和峰值移動隨時間的變化情況.同時將相鄰2個時間窗內的相位均值相減，分析其均值差異隨時間的變化情況.

3 結果與討論

3.1 不同運動階段之間的鎖相相位分布差異

如圖3、圖4所示，在運動的各個階段，鎖相相位呈現出了差異性的變化.獼猴A中，在M1腦區等待階段，鎖相的相位分布出現了雙峰值的特點，雙峰位置分別出現在相位0.81和5.90處；運動進入反應階段后，峰值位置向后移動，位于相位2.14處；在伸手階段，峰值在相位3.39處出現；在抓握階段，峰值相位出現在4.71處.

a、b、c、d.M1腦區；e、f、g、h.S1腦區；i、j、k、l.PRC腦區.圖3 獼猴A在不同運動階段的相位概率密度分布Fig.3 Phase probability distribution of monkey A at different stages

a、b、c、d.M1腦區；e、f、g、h.S1腦區；i、j、k、l.PRC腦區.圖4 獼猴B在不同運動階段的相位概率密度分布Fig.4 Phase probability density distribution of monkey B at different stages

在同一腦區,獼猴B中，4個階段呈現出相同的變化情況，但峰值位置有所不同，在等待階段，獼猴B的鎖相相位分布同樣具有2個峰值，分別位于相位0.57和相位5.59處；在反應階段，相位分布峰值向后移動，出現在了相位2.38處；在伸手階段，相位峰值在2.95附近；進入抓握階段后，相位峰值則位于4.15附近.峰值分布情況見表1.通過實驗發現，在連續運動的過程中，相位分布峰值的變化隨著運動的進行而不斷向大弧度值移動，在等待階段，3個腦區內都出現了雙峰值的情況，在獼猴B的S1腦區進入伸手階段，雖然也出現了雙峰值的情況，但是其2個峰值的位置分別為2.45和3.52，峰值較為接近，與等待階段不同.

表1 獼猴不同運動階段鎖相相位分布峰值

鎖相相位峰值的變化，也體現在了各階段相位均值的變化上，由圖5可以看出，由于鎖相相位在不同運動階段的分布變化，使得相位的均值呈現出了一定的變化規律，在2只獼猴上的3個腦區中都呈現出了一致的變化特點，從準備階段進入反應階段后，相位均值出現降低，隨后隨著運動的進行，到伸手階段時相位均值增加，到抓握階段時到達最大值.兩兩階段之間均進行過Wilcoxon秩和檢驗，檢驗結果P值均遠遠小于0.05，可以認為各階段之間具有明顯差異性.

圖5 不同階段鎖相相位均值變化Fig.5 Variation of mean phase of phase-locking at different stages

一個腦區鎖相相位的變化是該腦區所有通道鎖相相位變化的綜合體現，為了進一步分析腦區內鎖相相位分布變化是否具有一致性，對每一個腦區內的各個通道的相位均值進行了計算，并統計其均值變化與腦區整體變化一致的比例，如表2所示，可以看出所有腦區內的通道變化相同的比例均大于60%，在獼猴B的S1和PPC腦區，這一比例達到100%.

表2 均值變化與腦區整體變化相同的通道比例

3.2 事件發生前后100 ms窗口內的鎖相相位分布差異

不同運動階段由不同的事件進行劃分，在事件發生時，對鎖相相位變化的影響尚未知曉.為進一步分析運動過程中的鎖相相位變化，對各事件發生前后100 ms范圍內的鎖相相位分布進行了分析，如圖6、圖7所示，TargetOn、CenterRelease、TargetHit 3個事件前后各100 ms內分別對應等待階段和反應階段，反應階段和伸手階段，伸手階段和抓握階段.結果顯示，除獼猴B在CenterRelease事件前后，峰值無明顯變化，其他腦區在各個事件前后各100 ms內，相位的分布依然發生著變化，且隨著運動的進行，鎖相相位的峰值在向后移動，峰值情況見表3.所有事件前后的鎖相相位均經過置換檢驗，P值均遠小于0.05，可以認為在事件發生前后的100 ms里，鎖相相位分布發生了明顯變化.

圖6 事件發生前后100 ms內的獼猴A鎖相相位分布Fig.6 Plase-locking distribution of monkey A within 100 ms before and after the events

圖7 事件發生前后100 ms內的獼猴B鎖相相位分布Fig.7 Phase-locking distribution of monkay B within 100 ms before and after the events

表3 事件發生前后100 ms內鎖相相位分布峰值

在事件發生前后100 ms內的鎖相相位中，其均值同樣具有一些規律性的變化.如圖8所示，在2只獼猴的M1腦區，各個事件發生后的100 ms內的鎖相相位均值都要比事件發生前的100 ms內的鎖相相位均值要大；但在S1和PPC腦區內,TargetOn事件發生的前100 ms內的鎖相相位均值要比事件發生后的100 ms內的均值要大，即在等待階段將要結束的最后100 ms內的鎖相相位均值要比反應階段開始100 ms內的鎖相相位均值略高.

a、b、c.分別為獼猴A的M1、S1、PPC；d、e、f.分別為獼猴B的M1、S1、PPC.圖8 事件發生前后100 ms內的鎖相相位均值Fig.8 Mean value of the phase-locking phase in bins of 100 ms before and after the events

3.3 連續時間窗口下的鎖相相位分布變化

在事件發生前后100 ms內，鎖相相位發生了明顯變化，這種變化是否是隨著時間不斷進行仍然未知，基于此，本文在連續時間窗口中分析鎖相相位分布的變化情況，結果如圖9所示.圖9中TO、CR、TH分別是事件TargetOn、CenterRelease、TargetHit發生的時間.通過圖9可以發現，鎖相的相位分布變化是時刻都在發生的，但是在TargetOn事件發生后，即獼猴運動進入反應階段后，鎖相的相位分布開始有了更為明顯的峰值變化，在等待階段，雖然也有鎖相分布峰值的出現，但是其分布較為平緩，這與在不同運動階段討論鎖相相位分布的結果是一致的.

圖9 連續時間窗下鎖相相位概率密度分布變化Fig.9 Variation of phase-locking probability density distribution under continuous time bins

這種鎖相分布峰值的陡峭性，由圖9可知在CenterRelease發生后開始逐漸明顯，即在由反應階段到伸手階段時，鎖相的分布情況峰值更為明顯，在進入伸手階段后，鎖相相位分布峰值愈加明顯且峰值位置隨時間逐漸變化.

本文對每2個相鄰時間窗內的鎖相相位進行了均值相減，其隨時間變化的特征反應了鎖相相位隨時間變化的劇烈程度，如圖9所示.從圖10中可以看出在CenterRelease后突然增大的差異值，這代表了進入伸手階段后，鎖相相位變化速度的加劇.如表4所示，展示了2只獼猴各腦區出現最大差異值的時刻，該時刻表示了在本文所描述的獼猴伸展抓握實驗范式中，spike與LFP的鎖相相位變化最大的時刻.

圖10 連續相鄰時間窗的均值差變化Fig.10 Variation of mean difference between successive adjacent time bins

表4 連續相鄰時間窗的均值差最大值時間

4 結論

本文利用獼猴伸展抓握運動的實驗信號，在M1、S1、PPC腦區，對連續運動中的spike發放時刻和LFP delta節律的鎖相相位分布變化進行了研究.分別在運動的不同階段，關鍵事件發生前后100 ms內和連續時間窗3種時間尺度下，對鎖相相位的分布變化進行了分析.

獼猴在連續運動中的不同階段，spike與LFP delta節律的鎖相相位分布會發生顯著變化，在不同腦區的不同階段表現出了不同的相位分布特點和分布峰值位置.

為了進一步分析這種變化，對關鍵的事件TargetOn、CenterRelease、TargetHit發生前后各100 ms的時間窗內的鎖相相位分布進行了分析，可以發現事件前后2個相鄰100 ms內的鎖相相位分布是具有顯著差異性的，而且可以明顯觀測到鎖相相位的峰值隨時間推移逐漸增大.

為了細化分析這種峰值變化和時間的關系，本文使用連續時間窗對運動中各個時間點的鎖相相位分布進行了分析，實驗結果表明鎖相的相位分布隨時間不停變化，這種變化在CenterRelease前后最為明顯，相鄰窗口均值差異分析體現了相鄰時間窗之間的鎖相相位變化劇烈程度，可以看出在進入伸手階段后，這種變化最為劇烈，通過計算，得出了變化最大的時間點.

綜上所述，本文利用獼猴的伸展抓握運動范式，發現了峰電位與局部場電位delta節律鎖相相位分布的變化規律，所發現的鎖相相位的變化規律可以區分連續運動中的不同階段，為研究大腦神經信息編碼提供了一種新的思路.

致謝:本文實驗中所使用的獼猴伸展抓握運動范式數據由軍事醫學研究院軍事認知與腦科學研究所前沿交叉學科研究室提供，特別感謝該團隊對本項工作的支持.