神經元動作電位(鋒電位)的發放與局部場電位(LFP)節律之間的鎖相關系反映了重要的神經編碼信息,但目前的分析方法只能判斷有無鎖相性,而不能定量區別各種鎖相性的強弱。本文利用鋒電位觸發的疊加平均(STA)信號,以STA信號與原始LFP信號的功率百分比值φ為指標,來評價這種鎖相性的強弱。大鼠海馬CA1區的實驗數據以及仿真數據的驗證結果均表明,指標φ與鎖相性強弱之間存在單調關系,它能夠提供有效的界值來區分是否鎖相。由于計算該指標時無需事先濾波,它可以避免傳統的區間統計法中人為限定LFP頻率所造成的不利影響,為研究大腦神經信息編碼的機制提供一種新的定量方法。
引用本文： 吳蘊蘊, 封洲燕, 鄭曉靜. 定量分析神經元鋒電位與局部場電位鎖相關系的新方法. 生物醫學工程學雜志, 2014, 31(1): 172-180. doi: 10.7507/1001-5515.20140034 復制
引言
大腦神經元的主要電活動包括動作電位[即鋒電位(spike)]和突觸電位等，鋒電位是神經信息傳輸的基礎，而突觸電位等形成的局部場電位(local field potentials，LFP)在神經信息編碼和大腦認知等活動的研究中也呈現出越來越重要的作用[1-3]。
LFP中包含了各種節律波，頻率從低到高為：δ波(0.5~4 Hz)、θ波(4~8 Hz)、α波(8~12 Hz)、β波(12~30 Hz)和γ波(30~80 Hz)等。近年來的研究成果表明，鋒電位的發放與LFP節律之間的鎖相性(phase-locking)反映了重要的生理學信息。例如：大鼠嗅皮層錐體神經元的脈沖發放與LFP的γ波之間同步，存在緊密的相位耦合，γ振蕩波就像濾波器，能夠控制神經元的發放時間[4]；仿真研究表明神經元脈沖發放與LFP節律之間的鎖相關系可能編碼認知事件[5]；獼猴視皮層的γ節律可以編碼圖像的對比度信息[6]；甚至在人腦里也已發現，鋒電位發放時刻與LFP的θ節律之間的協調關系與記憶的形成密切相關[7]。因此，鋒電位與LFP節律之間鎖相性的研究具有重要意義。
判斷這兩者之間鎖相性的常用方法是區間統計法(簡稱Bin法)[8-9]，也就是先濾波提取LFP中的窄頻帶信號，如θ或者γ節律等，再把節律波的每個周期均分成若干個相位區間，統計神經元在每個區間上的發放概率，得到相位分布直方圖，最后利用統計檢驗來判斷神經元的發放與節律波之間是否存在顯著的鎖相性。
雖然Bin法的分布圖可以直觀顯示神經元的發放情況，但是，它有如下局限性：① 利用統計檢驗，只能判斷是否具有鎖相性，無法定量區分各種鎖相性的強弱。② 將LFP濾波成窄頻帶信號，這是一種人為的分割，有可能破壞LFP原本包含的重要信息。例如，大鼠清醒時θ節律的頻率為6~9 Hz，而睡眠和麻醉狀態時卻會變成2~5 Hz[9]。如果事先設定的濾波頻率范圍不合理，就可能損失信息。③ Bin法分析多種不同節律時需要重復統計工作，計算效率較低。另外，以相位分區的方式統計落在每個區間上的鋒電位個數時，區間大小的設定會直接影響統計檢驗的結果。
為了克服Bin法的上述局限性，本文利用鋒電位觸發的疊加平均(spike-triggered average，STA)法，直接從寬頻帶記錄信號中提取與鋒電位發放同步的場電位成分，并且設計定量指標來描述鋒電位與場電位節律波之間鎖相性的強弱。本文將以目前常見的θ和γ節律為例，用大鼠海馬CA1區記錄的實驗數據以及仿真數據來論證這種新方法評價鎖相性的有效性。
1 方法
1.1 動物實驗與信號采集
成年Sprague-Dawley大鼠(260～300 g，購自浙江省醫學科學院實驗動物中心)，腹腔注射烏拉坦(1.5 g/kg)麻醉之后，固定于大鼠腦立體定位儀上,切開頭部皮膚,除去部分顱骨。將記錄電極經大腦皮層插入海馬CA1區，定位是前囟后3.0 mm，旁開2.6 mm，由大腦皮層表面向下約2.5 mm。上下微調電極的位置，直到其中的某些通道記錄到明顯的鋒電位信號為止[10-11]。
記錄電極采用美國NeuroNexus Technologies公司生產的16通道微電極陣列。測量信號經過3600型16通道細胞外放大器(A-M Systems Inc. USA)放大100倍(頻率范圍為0.3 Hz ~ 5 kHz)，之后再使用PowerLab ML880型多通道信號采集系統(AD Instruments，Pty Ltd)采樣，頻率為20 kHz，A/D轉換分辨率為16位，并將獲得的原始采樣數據存入硬盤用于離線分析。這種原始信號稱寬頻帶信號，其中同時包含了LFP和多個神經元的鋒電位[10, 12]。16個通道可以同時記錄海馬CA1區從樹突至胞體不同結構層次上的神經電信號，本文選擇鋒電位幅值較大的胞體層通道用于檢測鋒電位，并選擇胞體層或者樹突層的通道用于分析LFP的θ和γ節律[9, 13]。
1.2 鋒電位的檢測和分類
利用PowerLab系統的LabChart軟件，將原始寬頻帶信號進行高通濾波，截止頻率設為500 Hz；然后，使用閾值法并手工輔助進行鋒電位檢測和分類，獲得各個鋒電位的負峰時間點作為發放時刻；并根據鋒電位的波形特征和發放特性[11, 13]，區分中間神經元和錐體神經元的鋒電位。其中，錐體神經元的鋒電位波形較寬，常以串的形式爆發式發放；而中間神經元的鋒電位波形較窄，通常均勻地發放。
1.3 鋒電位與場電位節律之間鎖相關系的分析
1.3.1 疊加平均分析法
鋒電位分類之后，以高通濾波信號上的鋒電位負峰時刻為基準，對應到所選擇通道(不一定為鋒電位所在通道)的原始寬頻帶信號上，向前和向后各截取所需長度的信號段(本文設為1 s)；然后求同類n個鋒電位觸發的疊加平均信號STA[14-16]。此STA信號中保留了與鋒電位發放同步的場電位成分。
為了與Bin法比較，進一步利用帶通濾波器對STA信號濾波，提取其中包含的所需頻段的場電位，如2~5 Hz的θ波和30～80 Hz的γ波等，分別稱為STA-θ信號和STA-γ信號。這樣，只要考察STA-θ和STA-γ信號中鋒電位所在的那個中心周期的波形，就可以求得鋒電位發放在場電位θ波和γ波節律上的相位均值。
當鋒電位個數n足夠大時，疊加平均可以消除LFP中包含的與鋒電位發放不同步的信號，保留并突顯出與鋒電位具有鎖相關系的LFP節律信號；并且，鎖相性越強，該節律信號在STA中的幅值就越大。但是，STA幅值的大小不僅與這種同步性的強弱有關，還與原始信號中節律波的原有幅值大小相關。因此，為了消除其影響，我們采用STA與原始LFP中感興趣節律的功率之比φ值，來評價鋒電位與LFP節律之間同步性的強弱，φ的定義如下。
根據帕塞瓦爾(Parserval)能量守恒定理，即信號f(t)的時域能量與頻域能量相等：
| $W=\int_{-\infty }^{\infty }{{{\left[ f\left( t \right) \right]}^{2}}dt}=\frac{1}{2\pi }\int_{-\infty }^{\infty }{{{\left| F\left( jw \right) \right|}^{2}}dw},$ |
而能量除以時間就得到功率。因此，信號的功率在時域或者頻域都可以計算。
由于STA-θ或者STA-γ信號中鋒電位所在的那個中心周期的信號較短，因此，在時域上計算其功率，即直接求數據的幅值平方和，再除以數據長度。而原始LFP信號較長，就在頻域上計算其包含的θ或者γ波的功率，方法是：利用LabChart軟件，用Welch周期圖法計算原始信號的功率譜密度。其中，FFT長度取65 536個數據點，對于20 kHz的采樣頻率，其頻譜分辨率約為0.305 Hz；窗函數采用Hanning窗[17]。然后求得θ波頻段或者γ波頻段的總功率。最后，得到STA與LFP的功率百分比φ，其計算公式為
| $\varphi =\frac{STA-\theta \left( 或STA-\gamma \right)中心周期波形幅值的平方和/該周期采樣點數}{LFP中包含的\theta \left( 或\gamma \right)頻段的功率}\times 100%,$ |
可以預見，φ值越大，表明鋒電位發放與LFP中節律波的鎖相性就越強。
1.3.2 相位區間分析法
為了檢驗上述STA法分析鋒電位發放與場電位節律之間鎖相關系的有效性，我們用相位區間法(即Bin法)[8, 18]分析同樣的實驗數據，并比較兩者的結果。在完成鋒電位檢測和分類的基礎上，Bin法的實現步驟如下：
① 確定節律波的各個周期
將原始寬頻帶信號進行帶通濾波，分別獲得θ波和γ波的頻段信號。然后，根據濾波信號上的所有過零點，將信號分段，各個分段要么包含波谷，要么包含波峰。相鄰兩個波峰之間的信號就作為節律波的一個周期，長度設為T。
② 計算鋒電位發放的相位分布
對于θ節律，若0.2 s≤T≤0.5 s，則為θ波；否則忽略此段信號以及該時間段內包含的鋒電位。將每個θ周期T(即0~360°)均分成18個相位區間(即Bin)，每個區間包含的角度為20°。統計已分類的同類各個鋒電位落在每個區間內的數目，計算發放概率，并作相位分布直方圖。
對于γ節律，與θ節律相同，測定30～80 Hz窄頻帶信號中每個γ波的周期。不過，γ波的出現是非連續的[13]，為了避免非γ節律信號段的影響，我們用閾值法選擇幅值足夠大的γ節律用于分析。閾值設定如下：用25 ms寬度的時間窗，將信號分段；計算各分段內所有數據點的幅值均方根，然后求這些均方根數值的均值及標準差，用均值加2倍標準差作為閾值[13]。如果某個周期波的谷點的幅值超過閾值，則選為γ波；否則忽略不計。本文分析的27段原長度為(4.7±4.4)min的信號中γ節律占24.4%±10.6%。最后，γ節律的相位分區和鋒電位發放的相位分布圖的作法都與θ節律的相同。
③ 確定平均相位角及其統計學意義
相位角數據是一種圓分布資料，要用圓分布統計法計算其均值和標準差，以便與STA法的結果進行比較。
假設某類n個鋒電位的相位角分別為α1，α2,…，αn令μ為平均角，則有如下計算公式[19]：
| $\left\{ \begin{array}{*{35}{l}} \bar{C}=\frac{1}{n}\sum\limits_{i=1}^{n}{\cos {{\alpha }_{i}}=\bar{r}\cos \mu } \\ \bar{S}=\frac{1}{n}\sum\limits_{i=1}^{n}{\sin {{\alpha }_{i}}=\bar{r}\sin \mu } \\ \bar{r}=\sqrt{{{{\bar{C}}}^{2}}+{{{\bar{S}}}^{2}}} \\ \end{array} \right.,$ |
其中
由上式可得，
| $\varpi =\pi \left\{ \begin{matrix} {{\tan }^{-1}}\left( \bar{S}/\bar{C} \right),若\bar{C}\ge 0 \\ {{\tan }^{-1}}\left( \bar{S}/\bar{C} \right)+\pi ,若\bar{C}＜0 \\ \end{matrix} \right.,$ |
設定μ取值在0～360°范圍，若計算得到的μ為負，則令μ=μ+360°。
αi的圓標準差(單位為弧度)是：
| $\sigma =\sqrt{-2\ln \bar{r}}\approx \sqrt{2\left( 1-\bar{r} \right)}$ |
若
| $P+\exp \left[ \sqrt{\left( 1+4n+4\left( {{n}^{2}}-{{\left( n\bar{r} \right)}^{2}} \right) \right)}-\left( 1+2n \right) \right]$ |
當P值大于所設定的檢驗水平(如0.05)時，就判斷為αi的分布不均勻[19]，也就是鋒電位的發放與所研究LFP節律的相關性具有統計學意義，本文簡稱鎖相，否則稱為不鎖相。
另外，本文在設計仿真數據時還使用了圓分布統計學中的正態分布，即von Mises分布[20]，其概率分布函數為
| $P\left( x;\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\mu },k \right)=\frac{1}{2\pi {{I}_{0}}\left( k \right)}{{e}^{k\cos \left( x-\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\mu } \right)}},$ |
式中x是自變量，I0是0階修正貝塞爾函數，參數和1/k分別相應于正態分布均值μ和方差σ2。
1.3.3 STA法與Bin法計算結果的相關性分析
本文分別使用上述STA法和Bin法計算了鋒電位發放相對于LFP中的θ和γ兩種節律的相位平均角，然后，利用圓相關來分析兩者之間的相關性。
圓相關系數R的計算公式為[20]：
| $R=\frac{\sum\limits_{i=1}^{n}{\sin \left( {{\alpha }_{i}}-{{\mu }_{a}} \right)\sin \left( {{\beta }_{i}}-{{\mu }_{\beta }} \right)}}{\sqrt{\sum\limits_{i=1}^{n}{{{\sin }^{2}}\left( {{\alpha }_{i}}-{{\mu }_{a}} \right){{\sin }^{2}}\left( {{\beta }_{i}}-{{\mu }_{\beta }} \right)}}},$ |
其中αi和βi分別為兩組樣本數據，μα和μβ分別是它們的平均角。假設兩組樣本之間不相關，則檢驗統計量為
| $t=\sqrt{f}\centerdot R$ |
它服從標準正態分布，其中
| $f=N\frac{\sum\limits_{i=1}^{n}{{{\sin }^{2}}\left( {{\alpha }_{i}}-{{\mu }_{a}} \right)\sum\limits_{i=1}^{n}{{{\sin }^{2}}\left( {{\beta }_{i}}-{{\mu }_{\beta }} \right)}}}{\sum\limits_{i=1}^{n}{{{\sin }^{2}}\left( {{\alpha }_{i}}-{{\mu }_{a}} \right){{\sin }^{2}}\left( {{\beta }_{i}}-{{\mu }_{\beta }} \right)}}$ |
2 結果
2.1 利用STA法提取與鋒電位發放同步的LFP成分
如圖 1(a)所示，在海馬CA1區記錄的兩路原始寬頻帶信號中同時包含了高頻鋒電位和低頻場電位，根據鋒電位波形和發放特性可知[11](見時間放大圖)，其中胞體層(pyramidal layer)通道含有明顯的錐體神經元，而樹突層(stratum radiatum)通道包含中間神經元鋒電位。經過鋒電位分類和疊加平均計算，這兩個神經元的STA信號如圖 1(b)所示，兩個STA中都有明顯的θ節律。分別將兩個STA的功率譜與原始信號功率譜進行比較(見圖 1c)，可見，STA的功率衰減很大，但它們分別在3.7 Hz和3.1 Hz處存在波峰，其它波峰則消失，因此，兩個神經元的發放可能都與θ節律具有鎖相性。
圖1
STA法
(a)在大鼠海馬CA1區胞體層和樹突層中，記錄到的兩路原始寬頻帶記錄信號中的小段示例及其時間放大信號，它們分別包含了中間神經元和錐體神經元的鋒電位；(b)分別從兩路原始信號中求得的中間神經元和錐體神經元的STA信號；(c)STA信號與原始信號功率譜的比較
Figure1. STA method(a) two short examples of wide frequency-band raw signals recorded respectively in pyramidal layer and stratum radiatum in rat hippocampal CA1 region,with the time expanded insets showing spikes of a pyramidal cell and a interneuron; (b) STA signals calculated from the two raw signals; (c) comparisons of power spectra between raw signals and STA signals
本文共計算了海馬CA1區記錄的12個錐體神經元和15個中間神經元的STA信號。以圖 2所示為例，在STA的帶通濾波信號STA-θ和STA-γ上，可以測得鄰近鋒電位的θ和γ節律的周期。它們的平均頻率分別為3.1 Hz和51.1 Hz(見表 1)。下面我們用1.3.1節定義的百分比值φ來評價各個神經元發放與場電位θ和γ節律鎖相性的強弱。
圖2
STA信號中的節律測量
將STA信號濾波后得到STA-
the measurements of
表1
CA1區錐體神經元和中間神經元的STA信號中所包含的兩種節律
Table1.
θ and γ oscillations included in the STA signals obtained from spikes of pyramidal cells and interneurons in CA1 region
2.2 鋒電位與LFP同步性強弱的評價指標
如圖 3所示是2個神經元的θ節律相位分布圖及其STA信號，為了便于考察節律性，每張圖都重復畫出了兩個周期的直方圖。其中圖 3(a)所示神經元的鎖相性較強，而圖 3(b)所示神經元的鎖相性則較弱。從相應的STA信號及其φ值和相位分布Rayleigh檢驗的P值可見，φ值較大的鋒電位發放具有顯著的鎖相性(P＜0.001)。
圖3
兩個神經元信號的鋒電位發放相位分布圖、Rayleigh檢驗結果、STA信號及φ值
(a)鎖相神經元;(b)不鎖相神經元。為了更好地說明周期性，本柱狀圖包括兩個周期的信號
Figure3. Spike θ phase distribution histograms of two neurons,Rayleigh test results and their STA signals with φ-values(a) a phase-locked neuron; (b) a non-phase-locked neuron. The histograms are repeated for two cycles for clarity
對于上一節所述共27個神經元的數據，分別計算節律波較明顯的近樹突層通道的STA-θ和STA-γ信號的φ值，結果如圖 4所示。圖中用雙坐標分別表示兩種節律的φ值，其均值和標準差如表 2所示。其中鎖相性的判斷由鋒電位發放分布圖的Rayleigh檢驗獲得。
圖4
神經元θ節律和γ節律的STA的φ值(n=27)
其中鎖相性的判斷由鋒電位發放分布數據的Rayleigh檢驗獲得，
Diamonds and squares represented phase-locked and non-phase-locked data of Rayleigh tests,respectively. Arrows indicate three data that their
表2
神經元的STA-θ信號和STA-γ信號的φ值(對于θ波，如果以φ=0.2%為界，那么，除了圖 4箭頭所指的1個數據以外，其余26個數據中，當φ＜0.2%時不鎖相，反之鎖相；對于γ波，如果以φ=0.4%為界，那么，除了箭頭所指的2個數據以外，其余25個數據中，當φ＜0.4%時不鎖相，反之鎖相。
φ值與Rayleigh檢驗結果不一致的3個數據點對應2組數據源(其中1組數據同時存在θ和γ節律的不一致)，它們的鋒電位個數n偏少，分別為124和190個，遠小于錐體神經元和中間神經元的平均個數1206和2010(見表 1)，本文下面的仿真研究將說明其中的原因。
雖然利用Rayleigh檢驗得到的P值可以判斷鋒電位發放的鎖相性是否具有統計學意義，但是，該方法無法提供鎖相性強弱的信息。 而φ值反映了經過疊加平均之后STA信號中仍然保留的某種節律波的百分比含量，其數值的大小可以定量反映鎖相性的強弱，本文下面設計仿真數據進一步研究此定量關系。
用正弦波y=sin t模擬LFP中的節律波，并假設鋒電位的相位分布符合圓正態分布，即von Mises分布VM(μ,k,n)。令鋒電位個數n=1 000；相位均值μ=250°；k=1/σ2，σ表示相位的標準差，即鋒電位發放的相位集中程度。 σ很大時，相位分布趨于均勻。 如圖 5(a)所示，仿真信號STA的φ值隨σ的增加呈現單調下降，根據Rayleigh檢驗結果，當φ＜0.16%(σ＞79.3°)時，相位分布不鎖相。該φ值與實驗結果接近。
圖5
用Von-mises分布仿真鋒電位發放結果
(a)STA法的
(a) changes of STA
為了考察鋒電位個數n對于φ值和Rayleigh檢驗結果的影響，令von Mises分布參數μ=250°，σ=70°，依次取n=10、100、500、1 000、2 000，分別產生50組隨機數用于仿真鋒電位的相位。各個n取值下隨機數據STA的φ值和Rayleigh檢驗結果如圖 5(b)所示，圖中下方標出了各組隨機數實際的μ和σ的平均值，隨機性導致它們與設定值之間存在差別。圖示仿真結果表明，當n≥500時，φ值趨于穩定；當n太小時，φ值的變動很大(即標準差很大)。其中，n=10時，有40組數據Rayleigh檢驗結果為不鎖相，但它們具有較大的φ值8.1% ±6.9%；而n=100時，也有6組為不鎖相數據的φ值達到1.3%±0.9%。這說明圖 4所示個別實驗數據φ值大卻不鎖相的原因就是鋒電位個數太少。
2.3 STA法與Bin法計算的平均相位角之間的相關性
如圖 6所示為STA法和Bin法計算的鋒電位相位平均角之間的相關性(已去除圖 4中Rayleigh檢驗不鎖相的數據點)。可見，相對于LFP中的θ和γ節律，兩種方法得到的相位角之間都具有很好的相關性，圓相關系數的t檢驗都有P＜0.001，因此，STA法計算得到的相位角可以表示鋒電位發放在場電位節律上的鎖相關系。
圖6
用STA法和Bin法計算的鋒電位相位平均角之間的相關性
(a)
(a)
圖 6中用三角形和圓形分別表示錐體和中間神經元的鋒電位發放，它們的相位角均值和標準差如表 3所示。由于圓分布的0°與360°等價(本文定義為波峰)，因此，圖中數據表明錐體神經元在樹突層LFP的θ節律波峰和γ節律波谷附近發放較多，而中間神經元則相反。
表3
用STA法與Bin法計算錐體神經元與中間神經元信號中，鋒電位發放在θ節律和γ節律中對應的相位(圓周平均±圓周標準差)
Table3.
Comparison of spike phases (circular mean ± standard deviation) in θ and γ oscillations obtained by Bin method and STA method for pyramidal cells and interneurons
3 討論
本文利用鋒電位觸發的疊加平均法來研究神經元發放與場電位節律之間的鎖相關系，并提出了描述鎖相性強弱的定量指標φ百分比值，這種新方法與常用的Bin法相比，具有如下優點：
首先，雖然本文為了驗證方法的有效性，采用了常見的θ和γ節律的頻率段，但是，可以看出該疊加平均法的重要特點之一是可以直接從STA信號考察與鋒電位具有鎖相關系的LFP節律，以免人為濾波造成的影響。而且，可以同時考察不同頻率的多種節律波，無需重復進行統計計算，可以提高效率。
其次，本文提出的φ指標與鋒電位發放的集中程度之間存在單調對應關系，可以定量表征鎖相性的強弱。而且，對于海馬CA1區θ節律、γ節律以及單頻率正弦仿真波的研究結果表明，當鋒電位個數n足夠大時，Rayleigh檢驗鎖相與不鎖相所對應的φ界值很接近，分別為0.4%、0.2%和0.16%，可用于判斷鎖相性。當n較小時，不足以充分衰減STA中的非鎖相信號，因此，φ值偏大。同樣，Bin法的相位分布圖也需要足夠大的n，否則Rayleigh統計檢驗也會失效。
另外，雖然已有研究人員利用疊加平均的STA信號來比較神經元發放與LFP節律之間的關系在不同行為狀態下的變化[14-16]，但是，尚未見利用定量指標來區別鎖相性強弱的文獻報道。傳統的Bin法[8, 18]利用統計學檢驗，只能判斷是否存在鎖相性，而無法區分各種鎖相性的強弱。本文利用STA信號計算的φ指標解決了這個問題。
最后，本文應用φ指標研究麻醉大鼠CA1區神經元鋒電位發放的結果表明，錐體神經元在樹突層LFP的θ節律波峰和γ節律波谷附近發放較多，而中間神經元則相反。已知樹突層與胞體層的節律波具有反相位關系[9]，因此，此研究結果與已有的報道一致[21-22]。
總之，本文利用疊加平均法研究神經元鋒電位發放與場電位振蕩節律之間的同步關系，所提出的定量指標能夠有效區分兩者之間鎖相性的強弱，從而為研究神經信息編碼等大腦活動的機制提供了一種新方法。
引言
大腦神經元的主要電活動包括動作電位[即鋒電位(spike)]和突觸電位等，鋒電位是神經信息傳輸的基礎，而突觸電位等形成的局部場電位(local field potentials，LFP)在神經信息編碼和大腦認知等活動的研究中也呈現出越來越重要的作用[1-3]。
LFP中包含了各種節律波，頻率從低到高為：δ波(0.5~4 Hz)、θ波(4~8 Hz)、α波(8~12 Hz)、β波(12~30 Hz)和γ波(30~80 Hz)等。近年來的研究成果表明，鋒電位的發放與LFP節律之間的鎖相性(phase-locking)反映了重要的生理學信息。例如：大鼠嗅皮層錐體神經元的脈沖發放與LFP的γ波之間同步，存在緊密的相位耦合，γ振蕩波就像濾波器，能夠控制神經元的發放時間[4]；仿真研究表明神經元脈沖發放與LFP節律之間的鎖相關系可能編碼認知事件[5]；獼猴視皮層的γ節律可以編碼圖像的對比度信息[6]；甚至在人腦里也已發現，鋒電位發放時刻與LFP的θ節律之間的協調關系與記憶的形成密切相關[7]。因此，鋒電位與LFP節律之間鎖相性的研究具有重要意義。
判斷這兩者之間鎖相性的常用方法是區間統計法(簡稱Bin法)[8-9]，也就是先濾波提取LFP中的窄頻帶信號，如θ或者γ節律等，再把節律波的每個周期均分成若干個相位區間，統計神經元在每個區間上的發放概率，得到相位分布直方圖，最后利用統計檢驗來判斷神經元的發放與節律波之間是否存在顯著的鎖相性。
雖然Bin法的分布圖可以直觀顯示神經元的發放情況，但是，它有如下局限性：① 利用統計檢驗，只能判斷是否具有鎖相性，無法定量區分各種鎖相性的強弱。② 將LFP濾波成窄頻帶信號，這是一種人為的分割，有可能破壞LFP原本包含的重要信息。例如，大鼠清醒時θ節律的頻率為6~9 Hz，而睡眠和麻醉狀態時卻會變成2~5 Hz[9]。如果事先設定的濾波頻率范圍不合理，就可能損失信息。③ Bin法分析多種不同節律時需要重復統計工作，計算效率較低。另外，以相位分區的方式統計落在每個區間上的鋒電位個數時，區間大小的設定會直接影響統計檢驗的結果。
為了克服Bin法的上述局限性，本文利用鋒電位觸發的疊加平均(spike-triggered average，STA)法，直接從寬頻帶記錄信號中提取與鋒電位發放同步的場電位成分，并且設計定量指標來描述鋒電位與場電位節律波之間鎖相性的強弱。本文將以目前常見的θ和γ節律為例，用大鼠海馬CA1區記錄的實驗數據以及仿真數據來論證這種新方法評價鎖相性的有效性。
1 方法
1.1 動物實驗與信號采集
成年Sprague-Dawley大鼠(260～300 g，購自浙江省醫學科學院實驗動物中心)，腹腔注射烏拉坦(1.5 g/kg)麻醉之后，固定于大鼠腦立體定位儀上,切開頭部皮膚,除去部分顱骨。將記錄電極經大腦皮層插入海馬CA1區，定位是前囟后3.0 mm，旁開2.6 mm，由大腦皮層表面向下約2.5 mm。上下微調電極的位置，直到其中的某些通道記錄到明顯的鋒電位信號為止[10-11]。
記錄電極采用美國NeuroNexus Technologies公司生產的16通道微電極陣列。測量信號經過3600型16通道細胞外放大器(A-M Systems Inc. USA)放大100倍(頻率范圍為0.3 Hz ~ 5 kHz)，之后再使用PowerLab ML880型多通道信號采集系統(AD Instruments，Pty Ltd)采樣，頻率為20 kHz，A/D轉換分辨率為16位，并將獲得的原始采樣數據存入硬盤用于離線分析。這種原始信號稱寬頻帶信號，其中同時包含了LFP和多個神經元的鋒電位[10, 12]。16個通道可以同時記錄海馬CA1區從樹突至胞體不同結構層次上的神經電信號，本文選擇鋒電位幅值較大的胞體層通道用于檢測鋒電位，并選擇胞體層或者樹突層的通道用于分析LFP的θ和γ節律[9, 13]。
1.2 鋒電位的檢測和分類
利用PowerLab系統的LabChart軟件，將原始寬頻帶信號進行高通濾波，截止頻率設為500 Hz；然后，使用閾值法并手工輔助進行鋒電位檢測和分類，獲得各個鋒電位的負峰時間點作為發放時刻；并根據鋒電位的波形特征和發放特性[11, 13]，區分中間神經元和錐體神經元的鋒電位。其中，錐體神經元的鋒電位波形較寬，常以串的形式爆發式發放；而中間神經元的鋒電位波形較窄，通常均勻地發放。
1.3 鋒電位與場電位節律之間鎖相關系的分析
1.3.1 疊加平均分析法
鋒電位分類之后，以高通濾波信號上的鋒電位負峰時刻為基準，對應到所選擇通道(不一定為鋒電位所在通道)的原始寬頻帶信號上，向前和向后各截取所需長度的信號段(本文設為1 s)；然后求同類n個鋒電位觸發的疊加平均信號STA[14-16]。此STA信號中保留了與鋒電位發放同步的場電位成分。
為了與Bin法比較，進一步利用帶通濾波器對STA信號濾波，提取其中包含的所需頻段的場電位，如2~5 Hz的θ波和30～80 Hz的γ波等，分別稱為STA-θ信號和STA-γ信號。這樣，只要考察STA-θ和STA-γ信號中鋒電位所在的那個中心周期的波形，就可以求得鋒電位發放在場電位θ波和γ波節律上的相位均值。
當鋒電位個數n足夠大時，疊加平均可以消除LFP中包含的與鋒電位發放不同步的信號，保留并突顯出與鋒電位具有鎖相關系的LFP節律信號；并且，鎖相性越強，該節律信號在STA中的幅值就越大。但是，STA幅值的大小不僅與這種同步性的強弱有關，還與原始信號中節律波的原有幅值大小相關。因此，為了消除其影響，我們采用STA與原始LFP中感興趣節律的功率之比φ值，來評價鋒電位與LFP節律之間同步性的強弱，φ的定義如下。
根據帕塞瓦爾(Parserval)能量守恒定理，即信號f(t)的時域能量與頻域能量相等：
| $W=\int_{-\infty }^{\infty }{{{\left[ f\left( t \right) \right]}^{2}}dt}=\frac{1}{2\pi }\int_{-\infty }^{\infty }{{{\left| F\left( jw \right) \right|}^{2}}dw},$ |
而能量除以時間就得到功率。因此，信號的功率在時域或者頻域都可以計算。
由于STA-θ或者STA-γ信號中鋒電位所在的那個中心周期的信號較短，因此，在時域上計算其功率，即直接求數據的幅值平方和，再除以數據長度。而原始LFP信號較長，就在頻域上計算其包含的θ或者γ波的功率，方法是：利用LabChart軟件，用Welch周期圖法計算原始信號的功率譜密度。其中，FFT長度取65 536個數據點，對于20 kHz的采樣頻率，其頻譜分辨率約為0.305 Hz；窗函數采用Hanning窗[17]。然后求得θ波頻段或者γ波頻段的總功率。最后，得到STA與LFP的功率百分比φ，其計算公式為
| $\varphi =\frac{STA-\theta \left( 或STA-\gamma \right)中心周期波形幅值的平方和/該周期采樣點數}{LFP中包含的\theta \left( 或\gamma \right)頻段的功率}\times 100%,$ |
可以預見，φ值越大，表明鋒電位發放與LFP中節律波的鎖相性就越強。
1.3.2 相位區間分析法
為了檢驗上述STA法分析鋒電位發放與場電位節律之間鎖相關系的有效性，我們用相位區間法(即Bin法)[8, 18]分析同樣的實驗數據，并比較兩者的結果。在完成鋒電位檢測和分類的基礎上，Bin法的實現步驟如下：
① 確定節律波的各個周期
將原始寬頻帶信號進行帶通濾波，分別獲得θ波和γ波的頻段信號。然后，根據濾波信號上的所有過零點，將信號分段，各個分段要么包含波谷，要么包含波峰。相鄰兩個波峰之間的信號就作為節律波的一個周期，長度設為T。
② 計算鋒電位發放的相位分布
對于θ節律，若0.2 s≤T≤0.5 s，則為θ波；否則忽略此段信號以及該時間段內包含的鋒電位。將每個θ周期T(即0~360°)均分成18個相位區間(即Bin)，每個區間包含的角度為20°。統計已分類的同類各個鋒電位落在每個區間內的數目，計算發放概率，并作相位分布直方圖。
對于γ節律，與θ節律相同，測定30～80 Hz窄頻帶信號中每個γ波的周期。不過，γ波的出現是非連續的[13]，為了避免非γ節律信號段的影響，我們用閾值法選擇幅值足夠大的γ節律用于分析。閾值設定如下：用25 ms寬度的時間窗，將信號分段；計算各分段內所有數據點的幅值均方根，然后求這些均方根數值的均值及標準差，用均值加2倍標準差作為閾值[13]。如果某個周期波的谷點的幅值超過閾值，則選為γ波；否則忽略不計。本文分析的27段原長度為(4.7±4.4)min的信號中γ節律占24.4%±10.6%。最后，γ節律的相位分區和鋒電位發放的相位分布圖的作法都與θ節律的相同。
③ 確定平均相位角及其統計學意義
相位角數據是一種圓分布資料，要用圓分布統計法計算其均值和標準差，以便與STA法的結果進行比較。
假設某類n個鋒電位的相位角分別為α1，α2,…，αn令μ為平均角，則有如下計算公式[19]：
| $\left\{ \begin{array}{*{35}{l}} \bar{C}=\frac{1}{n}\sum\limits_{i=1}^{n}{\cos {{\alpha }_{i}}=\bar{r}\cos \mu } \\ \bar{S}=\frac{1}{n}\sum\limits_{i=1}^{n}{\sin {{\alpha }_{i}}=\bar{r}\sin \mu } \\ \bar{r}=\sqrt{{{{\bar{C}}}^{2}}+{{{\bar{S}}}^{2}}} \\ \end{array} \right.,$ |
其中
由上式可得，
| $\varpi =\pi \left\{ \begin{matrix} {{\tan }^{-1}}\left( \bar{S}/\bar{C} \right),若\bar{C}\ge 0 \\ {{\tan }^{-1}}\left( \bar{S}/\bar{C} \right)+\pi ,若\bar{C}＜0 \\ \end{matrix} \right.,$ |
設定μ取值在0～360°范圍，若計算得到的μ為負，則令μ=μ+360°。
αi的圓標準差(單位為弧度)是：
| $\sigma =\sqrt{-2\ln \bar{r}}\approx \sqrt{2\left( 1-\bar{r} \right)}$ |
若
| $P+\exp \left[ \sqrt{\left( 1+4n+4\left( {{n}^{2}}-{{\left( n\bar{r} \right)}^{2}} \right) \right)}-\left( 1+2n \right) \right]$ |
當P值大于所設定的檢驗水平(如0.05)時，就判斷為αi的分布不均勻[19]，也就是鋒電位的發放與所研究LFP節律的相關性具有統計學意義，本文簡稱鎖相，否則稱為不鎖相。
另外，本文在設計仿真數據時還使用了圓分布統計學中的正態分布，即von Mises分布[20]，其概率分布函數為
| $P\left( x;\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\mu },k \right)=\frac{1}{2\pi {{I}_{0}}\left( k \right)}{{e}^{k\cos \left( x-\overset{\lower0.5em\hbox{$\smash{\scriptscriptstyle\frown}$}}{\mu } \right)}},$ |
式中x是自變量，I0是0階修正貝塞爾函數，參數和1/k分別相應于正態分布均值μ和方差σ2。
1.3.3 STA法與Bin法計算結果的相關性分析
本文分別使用上述STA法和Bin法計算了鋒電位發放相對于LFP中的θ和γ兩種節律的相位平均角，然后，利用圓相關來分析兩者之間的相關性。
圓相關系數R的計算公式為[20]：
| $R=\frac{\sum\limits_{i=1}^{n}{\sin \left( {{\alpha }_{i}}-{{\mu }_{a}} \right)\sin \left( {{\beta }_{i}}-{{\mu }_{\beta }} \right)}}{\sqrt{\sum\limits_{i=1}^{n}{{{\sin }^{2}}\left( {{\alpha }_{i}}-{{\mu }_{a}} \right){{\sin }^{2}}\left( {{\beta }_{i}}-{{\mu }_{\beta }} \right)}}},$ |
其中αi和βi分別為兩組樣本數據，μα和μβ分別是它們的平均角。假設兩組樣本之間不相關，則檢驗統計量為
| $t=\sqrt{f}\centerdot R$ |
它服從標準正態分布，其中
| $f=N\frac{\sum\limits_{i=1}^{n}{{{\sin }^{2}}\left( {{\alpha }_{i}}-{{\mu }_{a}} \right)\sum\limits_{i=1}^{n}{{{\sin }^{2}}\left( {{\beta }_{i}}-{{\mu }_{\beta }} \right)}}}{\sum\limits_{i=1}^{n}{{{\sin }^{2}}\left( {{\alpha }_{i}}-{{\mu }_{a}} \right){{\sin }^{2}}\left( {{\beta }_{i}}-{{\mu }_{\beta }} \right)}}$ |
2 結果
2.1 利用STA法提取與鋒電位發放同步的LFP成分
如圖 1(a)所示，在海馬CA1區記錄的兩路原始寬頻帶信號中同時包含了高頻鋒電位和低頻場電位，根據鋒電位波形和發放特性可知[11](見時間放大圖)，其中胞體層(pyramidal layer)通道含有明顯的錐體神經元，而樹突層(stratum radiatum)通道包含中間神經元鋒電位。經過鋒電位分類和疊加平均計算，這兩個神經元的STA信號如圖 1(b)所示，兩個STA中都有明顯的θ節律。分別將兩個STA的功率譜與原始信號功率譜進行比較(見圖 1c)，可見，STA的功率衰減很大，但它們分別在3.7 Hz和3.1 Hz處存在波峰，其它波峰則消失，因此，兩個神經元的發放可能都與θ節律具有鎖相性。
圖1
STA法
(a)在大鼠海馬CA1區胞體層和樹突層中，記錄到的兩路原始寬頻帶記錄信號中的小段示例及其時間放大信號，它們分別包含了中間神經元和錐體神經元的鋒電位；(b)分別從兩路原始信號中求得的中間神經元和錐體神經元的STA信號；(c)STA信號與原始信號功率譜的比較
Figure1. STA method(a) two short examples of wide frequency-band raw signals recorded respectively in pyramidal layer and stratum radiatum in rat hippocampal CA1 region,with the time expanded insets showing spikes of a pyramidal cell and a interneuron; (b) STA signals calculated from the two raw signals; (c) comparisons of power spectra between raw signals and STA signals
本文共計算了海馬CA1區記錄的12個錐體神經元和15個中間神經元的STA信號。以圖 2所示為例，在STA的帶通濾波信號STA-θ和STA-γ上，可以測得鄰近鋒電位的θ和γ節律的周期。它們的平均頻率分別為3.1 Hz和51.1 Hz(見表 1)。下面我們用1.3.1節定義的百分比值φ來評價各個神經元發放與場電位θ和γ節律鎖相性的強弱。
圖2
STA信號中的節律測量
將STA信號濾波后得到STA-
the measurements of
表1
CA1區錐體神經元和中間神經元的STA信號中所包含的兩種節律
Table1.
θ and γ oscillations included in the STA signals obtained from spikes of pyramidal cells and interneurons in CA1 region
2.2 鋒電位與LFP同步性強弱的評價指標
如圖 3所示是2個神經元的θ節律相位分布圖及其STA信號，為了便于考察節律性，每張圖都重復畫出了兩個周期的直方圖。其中圖 3(a)所示神經元的鎖相性較強，而圖 3(b)所示神經元的鎖相性則較弱。從相應的STA信號及其φ值和相位分布Rayleigh檢驗的P值可見，φ值較大的鋒電位發放具有顯著的鎖相性(P＜0.001)。
圖3
兩個神經元信號的鋒電位發放相位分布圖、Rayleigh檢驗結果、STA信號及φ值
(a)鎖相神經元;(b)不鎖相神經元。為了更好地說明周期性，本柱狀圖包括兩個周期的信號
Figure3. Spike θ phase distribution histograms of two neurons,Rayleigh test results and their STA signals with φ-values(a) a phase-locked neuron; (b) a non-phase-locked neuron. The histograms are repeated for two cycles for clarity
對于上一節所述共27個神經元的數據，分別計算節律波較明顯的近樹突層通道的STA-θ和STA-γ信號的φ值，結果如圖 4所示。圖中用雙坐標分別表示兩種節律的φ值，其均值和標準差如表 2所示。其中鎖相性的判斷由鋒電位發放分布圖的Rayleigh檢驗獲得。
圖4
神經元θ節律和γ節律的STA的φ值(n=27)
其中鎖相性的判斷由鋒電位發放分布數據的Rayleigh檢驗獲得，
Diamonds and squares represented phase-locked and non-phase-locked data of Rayleigh tests,respectively. Arrows indicate three data that their
表2
神經元的STA-θ信號和STA-γ信號的φ值(對于θ波，如果以φ=0.2%為界，那么，除了圖 4箭頭所指的1個數據以外，其余26個數據中，當φ＜0.2%時不鎖相，反之鎖相；對于γ波，如果以φ=0.4%為界，那么，除了箭頭所指的2個數據以外，其余25個數據中，當φ＜0.4%時不鎖相，反之鎖相。
φ值與Rayleigh檢驗結果不一致的3個數據點對應2組數據源(其中1組數據同時存在θ和γ節律的不一致)，它們的鋒電位個數n偏少，分別為124和190個，遠小于錐體神經元和中間神經元的平均個數1206和2010(見表 1)，本文下面的仿真研究將說明其中的原因。
雖然利用Rayleigh檢驗得到的P值可以判斷鋒電位發放的鎖相性是否具有統計學意義，但是，該方法無法提供鎖相性強弱的信息。 而φ值反映了經過疊加平均之后STA信號中仍然保留的某種節律波的百分比含量，其數值的大小可以定量反映鎖相性的強弱，本文下面設計仿真數據進一步研究此定量關系。
用正弦波y=sin t模擬LFP中的節律波，并假設鋒電位的相位分布符合圓正態分布，即von Mises分布VM(μ,k,n)。令鋒電位個數n=1 000；相位均值μ=250°；k=1/σ2，σ表示相位的標準差，即鋒電位發放的相位集中程度。 σ很大時，相位分布趨于均勻。 如圖 5(a)所示，仿真信號STA的φ值隨σ的增加呈現單調下降，根據Rayleigh檢驗結果，當φ＜0.16%(σ＞79.3°)時，相位分布不鎖相。該φ值與實驗結果接近。
圖5
用Von-mises分布仿真鋒電位發放結果
(a)STA法的
(a) changes of STA
為了考察鋒電位個數n對于φ值和Rayleigh檢驗結果的影響，令von Mises分布參數μ=250°，σ=70°，依次取n=10、100、500、1 000、2 000，分別產生50組隨機數用于仿真鋒電位的相位。各個n取值下隨機數據STA的φ值和Rayleigh檢驗結果如圖 5(b)所示，圖中下方標出了各組隨機數實際的μ和σ的平均值，隨機性導致它們與設定值之間存在差別。圖示仿真結果表明，當n≥500時，φ值趨于穩定；當n太小時，φ值的變動很大(即標準差很大)。其中，n=10時，有40組數據Rayleigh檢驗結果為不鎖相，但它們具有較大的φ值8.1% ±6.9%；而n=100時，也有6組為不鎖相數據的φ值達到1.3%±0.9%。這說明圖 4所示個別實驗數據φ值大卻不鎖相的原因就是鋒電位個數太少。
2.3 STA法與Bin法計算的平均相位角之間的相關性
如圖 6所示為STA法和Bin法計算的鋒電位相位平均角之間的相關性(已去除圖 4中Rayleigh檢驗不鎖相的數據點)。可見，相對于LFP中的θ和γ節律，兩種方法得到的相位角之間都具有很好的相關性，圓相關系數的t檢驗都有P＜0.001，因此，STA法計算得到的相位角可以表示鋒電位發放在場電位節律上的鎖相關系。
圖6
用STA法和Bin法計算的鋒電位相位平均角之間的相關性
(a)
(a)
圖 6中用三角形和圓形分別表示錐體和中間神經元的鋒電位發放，它們的相位角均值和標準差如表 3所示。由于圓分布的0°與360°等價(本文定義為波峰)，因此，圖中數據表明錐體神經元在樹突層LFP的θ節律波峰和γ節律波谷附近發放較多，而中間神經元則相反。
表3
用STA法與Bin法計算錐體神經元與中間神經元信號中，鋒電位發放在θ節律和γ節律中對應的相位(圓周平均±圓周標準差)
Table3.
Comparison of spike phases (circular mean ± standard deviation) in θ and γ oscillations obtained by Bin method and STA method for pyramidal cells and interneurons
3 討論
本文利用鋒電位觸發的疊加平均法來研究神經元發放與場電位節律之間的鎖相關系，并提出了描述鎖相性強弱的定量指標φ百分比值，這種新方法與常用的Bin法相比，具有如下優點：
首先，雖然本文為了驗證方法的有效性，采用了常見的θ和γ節律的頻率段，但是，可以看出該疊加平均法的重要特點之一是可以直接從STA信號考察與鋒電位具有鎖相關系的LFP節律，以免人為濾波造成的影響。而且，可以同時考察不同頻率的多種節律波，無需重復進行統計計算，可以提高效率。
其次，本文提出的φ指標與鋒電位發放的集中程度之間存在單調對應關系，可以定量表征鎖相性的強弱。而且，對于海馬CA1區θ節律、γ節律以及單頻率正弦仿真波的研究結果表明，當鋒電位個數n足夠大時，Rayleigh檢驗鎖相與不鎖相所對應的φ界值很接近，分別為0.4%、0.2%和0.16%，可用于判斷鎖相性。當n較小時，不足以充分衰減STA中的非鎖相信號，因此，φ值偏大。同樣，Bin法的相位分布圖也需要足夠大的n，否則Rayleigh統計檢驗也會失效。
另外，雖然已有研究人員利用疊加平均的STA信號來比較神經元發放與LFP節律之間的關系在不同行為狀態下的變化[14-16]，但是，尚未見利用定量指標來區別鎖相性強弱的文獻報道。傳統的Bin法[8, 18]利用統計學檢驗，只能判斷是否存在鎖相性，而無法區分各種鎖相性的強弱。本文利用STA信號計算的φ指標解決了這個問題。
最后，本文應用φ指標研究麻醉大鼠CA1區神經元鋒電位發放的結果表明，錐體神經元在樹突層LFP的θ節律波峰和γ節律波谷附近發放較多，而中間神經元則相反。已知樹突層與胞體層的節律波具有反相位關系[9]，因此，此研究結果與已有的報道一致[21-22]。
總之，本文利用疊加平均法研究神經元鋒電位發放與場電位振蕩節律之間的同步關系，所提出的定量指標能夠有效區分兩者之間鎖相性的強弱，從而為研究神經信息編碼等大腦活動的機制提供了一種新方法。
