基于局部场电位固有模态分量的响应调谐特性研究

万 红 张晓娜 刘新玉 李晓燕

(郑州大学电气工程学院，郑州 450001)

引言

调谐特性是初级视觉皮层(primary visual cortex，又称V1 区)神经元感受野的主要特性之一［1］。传统的对V1 区神经元对空间频率响应特性的研究，大都基于神经元锋电位(spike)发放率来获得神经元对刺激响应的调谐曲线。然而在微电极阵列胞外采集技术中，spike 受到噪声的干扰经常会出现漏检和误检，即假阴性和假阳性问题，使后续分析受到严重影响。因此，与spike 活动同时记录到的局部场电位信号，逐渐成为目前视觉信息处理机制研究的热点之一［2］。

局部场电位(local field potential，LFP)是记录电极尖端附近局部区域的神经元兴奋性和抑制性突触后电位的总和。近年来，LFP 对于视觉刺激特征的调谐作用受到广泛关注。Gawne 等指出，局部场电位信号与锋电位发放率对于视觉刺激有着相似的调制作用，当锋电位发放率达到最大时，局部场电位的能量同时也达到了峰值［3］;Ince 等通过对朝向、空间频率、时间频率等3 个视觉刺激特征的研究，表明V1 区LFP 频带能量能够对上述特征进行有效编码［4］;Belitski 等也发现LFP 信号不同频带编码了不同的视觉刺激特征［5］。诸多研究表明LFP对不同视觉刺激特征具有调谐功能，然而对于LFP中与刺激的最佳响应范围，研究者们没有达成共识。

Henrie 等通过对麻醉猴的研究发现，LFP 在刺激前后功率谱会发生明显的变化，并且不同频带的能量变化不同，能量变化最明显频带主要集中在25～90 Hz［6］;Kavser 等以清醒猫为对象，研究了LFP不同频段对时间频率、空间频率和朝向的选择性，结果表明具有调谐特征的频带，主要集中在8 ～23 Hz 和39 ～109 Hz［7］;Philipp 等通过记录清醒猕猴主视区的LFP 信号，研究了神经元感受野的特征选择性，发现LFP 最佳调谐范围主要集中在30 ～90 Hz［8］。Magri 等通过记录麻醉猕猴初级视觉皮层的信号，将LFP 划分为表征自然视觉刺激的多个频带［9］。就目前研究结果来看，由于动物类型差异和LFP 的非平稳性，使得不同的实验获得的最佳响应范围并不一致。因此，获得LFP 最佳响应频带范围、提取出神经元响应特征，还需要深入研究。

局部场电位是典型的非线性、非平稳信号，经典的时频方法处理效果并不理想。近些年来，具有自适应特性的希尔伯特黄变换(Hilbert-Huang transform，HHT)正逐渐得到研究者的青睐。HHT 没有先验的基函数，而且分解过程完全基于数据本身［10］。Liang 等的研究表明，HHT 是神经电信号分析的强有力工具［11-13］。本研究基于HHT 算法，首先将LFP 信号进行经验模态分解，将其第二阶固有模态分量进行Hilbert 变换并求取其能量，利用该能量特征研究LFP 对不同空间频率刺激的调谐特性，并与同一电极记录到的多神经元锋电位(multi-unit activity，MUA)和小波分解提取的Gamma 频带进行了对比。

1 材料和方法

1.1 材料

实验动物选用年龄12 周左右、体重200 ～250 g的6 只成年Long Evans(LE)大鼠，雌雄不拘。术前腹腔注射10%水合氯醛麻醉，实施开颅手术。参考大鼠脑立体定位图谱，确定初级视觉皮层的中心位置，用颅钻进行开颅手术。用显微操纵器(MX7600，Siskiyou Inc.)植入2 ×8 微电极阵列(Microprobe，铂铱合金，电极间距250 μm，电极直径125 μm，尖端直径15 μm，电极阻抗在0.5 ～1.0 MΩ 之间)，并用数据采集系统(CerebusTM，Blockrock Inc.)记录大鼠V1 区神经元诱发放电信号，包括检测到的LFP 和MUA 信号。MUA 将单个通道的spike 作为一个整体，是电极尖端附近多个神经元胞体的活动［14］。视觉刺激采用不同空间频率的全屏运动光栅，空间频率分别取0.01、0.02、0.05、0.1、0.2、0.4 Cpd，时间频率为2 Hz，刺激1 s，间隔1 s。每组刺激各空间频率光栅随机出现，每次实验重复20 组。

1.2 方法

1.2.1 希尔伯特黄变换

HHT 包括经验模态分解(empirical mode decomposition，EMD)和Hilbert 变换两个部分。首先采用EMD 方法将LFP 信号分解为若干固有模态分量(intrinsic mode function，IMF)，然后对每个IMF 进行Hilbert 变换，求出其瞬时频率和幅值，构成被分析信号的时间-频率-幅值(能量)的三维谱图，即Hilbert 谱。EMD 分解的IMF 须满足以下两个条件:一是在整个信号长度上，极值点和过零点数必须相等或至多相差一个;二是在任何一点，由局部极大值点形成的包络线和局部极小值点形成的包络线的均值为零。

LFP 的HHT 变换过程如下:

首先，对其中一个通道LFP 信号x 进行EMD 分解，得到从高频到低频的n 个IMF 分量c1，c2，…，cn和余项rn，即

然后对每个IMF 分量进行Hilbert 变换

当y(t)与c(t)形成一个复共轭时，就可构造解析信号z(t)

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

从而得到IMF 分量的瞬时频率ωi(t)

式(4)和式(5)表明信号的幅值和瞬时频率都是时间的函数，则原始信号可表示为

取其实部，即原始信号的Hilbert 谱，

式中，Re 表示实部，这里省略了残余函数rn，因为它是一个单调函数或者常量，代表着长期振动。

1.2.2 响应特征提取

瞬时能量(instantaneous energy，IE)提供了一种能量随时间变化的信息，由Hilbert 谱可以得到

式中，H 为Hilbert 变换，ω 表示目标频带，t 代表时间。通过式(8)可以计算单一模态的瞬时能量。

1.2.3 调谐指数

调谐指数度量了信号特征对刺激参数的响应程度。为了定量的表示信号的调谐性能，分别计算了MUA 和LFP 在不同空间频率刺激下的调谐指数。MUA 的调谐指数dMUA为［8］

式中，Pmax表示MUA 发放率调谐曲线的最大值，Pmin表示MUA 发放率调谐曲线的最小值。同理，LFP 的调谐指数dLFP为［8］

式中，Emax表示LFP 能量调谐曲线的最大值，Emin表示LFP 能量调谐曲线的最小值。

为了定量衡量不同方法获得的神经元响应特征的一致程度，引入了神经元响应一致率指标，定义为:设总神经元个数为M，其中两种信号对最佳刺激空间频率响应一致的神经元个数为N，那么最佳空间频率刺激响应的一致率λ 定义为

1.2.4 小波分解提取LFP 的Gamma 频带

小波分解结果依赖于所选择的小波基。基于LFP 的特点，选择具有紧支正交基，满足精确信号重建条件的“db5”作为小波基函数［15］。小波分解表示为

小波重构

式中，小波基函数为

式中，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 进行小波重构，从而得到Gamma 频带(30 ～90 Hz)信号。

2 结果

2.1 各阶固有模态分量特征分析

首先对记录到的16 通道原始数据进行低通滤波，通带截止频率200 Hz，阶数200，利用零均值处理去除基线漂移，记为LFP。预处理前后的LFP 信号如图1 中(a)和(b)所示，其中刺激时间为1 s。将除余量之外的各阶IMF 信号合成求取Hilbert 谱，如图1(c)所示。由图可知，LFP 在刺激前后能量发生明显变化(t-test，P ＜0.05)。

为了进一步提取LFP 的响应频带范围，对LFP进行了EMD 分解，得到7 个固有模态分量和一个余量，如图1(d)所示。从分解结果来看，每一阶分量都具有不同振幅和频率，第一阶模态分量的频率最高，随着模态阶数的增加其所含的周期变长，直到最后一阶逼近一条线性的时间函数，而且前几个分量包含了信号的主要特征，在刺激后能量有不同程度的增加。

2.2 IMF2 的调谐性能

调谐特性的传统研究方法，大都基于spike 发放率，而且已得到了大量文献的证实［16-17］。因此选用MUA 发放率为基准，来统计分析各阶固有模态分量是否具有调谐特性。表1 给出了基于各阶固有模态分量，对不同空间频率的调谐指数平均值及与MUA 的一致率。从表中可以发现，LFP 的第二阶固有模态分量(IMF2)的调谐指数最大，其与MUA 最佳空间频率调谐的一致率为68.75%，显著高于其它固有模态分量。这一结果表明，LFP 携带的刺激相关信息主要集中在IMF2。

以IMF2 作为特征频带，以式(8)定义的计算单一模态的瞬时能量为能量特征，分析了大鼠V1 区空间频率响应特性，并将其与MUA 发放率进行了对比，结果如图2 所示。

图1 光栅刺激前后LFP 响应特性。(a)原始LFP 信号(刺激时间1 s);(b)预处理后LFP 信号;(c)LFP 的Hilbert谱;(d)LFP 经验模态分解后的各阶固有模态分量(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 固有模态分量调谐指数统计表Tab．1 The statistical table of tuning index based on the intrinsic mode functions of LFP

图2(a)表示20 次重复刺激IMF2 和MUA 发放率的对比结果。由图可知，IMF2 幅值和MUA 的发放率在刺激前后都发生了显著变化，两者具有一致性。图2(b)是20 次重复刺激IMF2 在灰屏和光栅刺激能量的对比图，虚线表示灰屏时的能量变化，实线表示光栅刺激时的能量变化，阴影部分为标准差，显著水平为0.05。由图可知光栅刺激下IMF2能量明显增加，进一步表明IMF2 携带了与刺激相关的信息。图2(c)为IMF2 和MUA 在不同空间频率刺激下的标准化调谐曲线(z-score 标准化，Matlab)，二者对不同空间频率刺激都有不同的响应强度，而且调谐曲线具有较好的一致性，同时在空间频率为0.02 Cpd 时响应最强烈。

2.3 IMF2 功率谱估计

为确定IMF2 的频带范围，对其进行了功率谱估计，并与小波分解提取的30 ～90 Hz 的Gamma 频带的功率谱进行了对比，如图3 所示。图中选用Welch 方法分别对IMF2 和Gamma 进行功率谱估计，Welch 估计的关键是选取合适的窗函数和窗长，选取主瓣较窄的Hamming 窗［15］。经对比研究发现:选用Hamming 窗分辨率20 Hz，即对应窗长N =200，能够准确提取频带能量。由图可知IMF2 的频带范围主要集中在35 ～80 Hz 内，比Gamma 频带范围更窄。

2.4 IMF2 与MUA、小波分解的调谐性能对比

图2 IMF2 与MUA 空间频率调谐特性对比。(a)20 次重复刺激下IMF2 和MUA 发放率;(b)灰屏和光栅刺激下IMF2 能量变化;(c)IMF2 与MUA 调谐曲线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

为进一步验证该方法在提取特征频带时的优越性，利用从6 只LE 大鼠记录到的共63 通道神经元响应信号，计算了IMF2 与MUA、小波分解提取的Gamma 频带对空间频率刺激响应具有一致性的通道调谐指数，结果如图4 所示。图4(a)为IMF2 和MUA 的调谐指数对比图。图中横坐标和纵坐标分别为IMF2 和MUA 的调谐指数，十字表示平均值。由图可知，IMF2 的平均调谐指数为0.795 1，MUA为0.631 3，这一结果表明IMF2 调谐性能要强于MUA。图4(b)是IMF2 与小波分解提取的Gamma频带的调谐指数对比图。图中横坐标和纵坐标分别为IMF2 和小波分解提取的Gamma 频带的调谐指数，十字表示平均值。由图可知，前者平均调谐指数为0.795 1，后者为0.664 6。由此可见，基于HHT方法提取的LFP 响应频带较小波分解方法更为有效。

3 讨论

LFP 活动反映了大脑皮层局部区域神经元集群的响应特性，可以从更精细的时间尺度和空间尺度上解释认知的过程，在视觉信息处理机制研究中得到了广泛应用。近年来研究发现，LFP 中与刺激相关的响应频带主要集中在Gamma 带，然而由于LFP的非平稳性，使得不同文献报道的响应频带范围并不一致。因此，确定LFP 的响应范围、提取LFP 的响应特征，还需要进一步深入研究。

图4 不同方法得到的调谐指数对比图。(a)IMF2 与MUA 调谐指数对比图;(b)IMF2 与小波分解的Gamma 频带调谐指数对比图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 的响应频带范围。由于HHT 变换引入了EMD 分解，使得经EMD 分解得到的各阶固有模态分量完全基于信号本身的特性，无需先验基函数，可以有效反应信号的固有特征，其分解的客观性、自适应性凸显了该方法在非平稳信号处理中的优势，具有广阔的应用前景。由实验结果可以看出，HHT 能够有效分解出LFP 的响应频带(如图1)，弥补了传统小波算法需要选取合适小波基的不足，拓展了LFP 处理的新思路，为LFP 信号的响应特征提取及后续研究奠定了基础。

本研究得出的结论与现有的研究结果比较吻合。目前的研究表明，LFP 的Gamma 频带(30 ～90 Hz)能量能够表征视觉刺激特征，本研究表明LFP经过EMD 分解得到的IMF2 的频带范围主要集中在35 ～70 Hz，这个范围的主要频带与Gamma 带重叠，而且较Gamma 频带范围更窄、更精确(如图3)，并且其调谐指数高于Gamma 频带(如图4(b))。这些结果表明基于HHT 方法提取LFP 响应信号更为有效。激的调谐特性，并与同一电极记录到的MUA 和小波分解提取的Gamma 频带进行了对比。

通过对记录的63 个通道神经元统计分析，LFP的IMF2 携带了刺激信息，与其他模态分量相比，IMF2 对刺激光栅空间频率的调谐特性最强，与MUA 空间频率响应的一致率达68.75%，且其调谐指数的均值(0.795 1)高于 MUA (0.631 3)和Gamma 频带(0.664 6)。表明基于LFP 第二固有模态分量来研究神经元感受野的最佳空间频率调谐特性，比MUA 和小波提取的Gamma 频带更有优势，为进一步从LFP 分析视觉皮层信息处理机制，提供了另一种可行的有效方法。

4 结论

本研究基于HHT 算法，首先将LFP 信号进行经验模态分解，然后将其第二阶固有模态分量作为特征频段，进行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］. 仪器仪表学报，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］ 尚志刚，冯平艳，刘新玉等. 局部场电位γ 频带能量对朝向调谐特性研究［J］. 郑州大学学报(工学版)，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.