基于前额联络皮层脑电信号的海洛因戒断大鼠个体觅药行为状态识别*

黄 磊， 潘群皖， 胡慧娴， 徐晓燕

((1. 皖南医学院 医学工程学教研室, 芜湖 241002; 2. 皖南医学院 生理教研室， 芜湖 241002)

药物成瘾是以强迫性觅药和复吸为特征性的慢性脑病[1]。研究表明，大脑前额叶联络皮层(frontal association cortex, FrA)中枢奖赏系统通路的正性强化作用，是导致大鼠海洛因成瘾的重要因素之一。Doherty等[2]研究显示，大鼠海洛因自给药戒断期，内侧前额皮层的边缘前区(prelimbic area, PL)和下缘区(infralimbic area, IL) Fos阳性神经元显著增加，这些神经元的活性与大鼠海洛因觅药行为正相关。场景暴露可激活这些Fos阳性神经元集群，其神经电活动编码海洛因奖赏效应与自给药场景联合觅药恢复过程[3]。

条件性位置偏爱(conditioned place preference，CPP)实验是测定动物药物成瘾的一种常用方法。使用FrA区电极埋藏加无线遥测技术，本课题组对大鼠阿片类药物位置偏爱状态下，FrA区神经元同步遥测脑电活动(局部场电位，local field potentials)进行分析，发现吗啡诱导CPP组大鼠黑-白箱(拌药箱)穿梭潜伏期及穿梭过程中，大鼠PL区小波提取的场电位各频段平均波幅均减小[4]，在使用海洛因(二乙酰吗啡)诱导CPP组大鼠实验中，也发现在黑-白箱穿梭潜伏期，双侧PL区场电位均呈现θ波百分比上升，左侧还表现为γ3波百分比增加，出现所谓的PL区θ波振荡(θ oscillation)的增强，尤其是左侧可能形成的θ～γ3同步振荡，呈现θ～γ波节律交错频率耦合(cross-frequency coupling，CFC)[5-7]等变化规律，同时发现表征脑电序列复杂度的熵值有变化[4-7]。这些均表明大鼠在黑-白箱穿梭觅药中，存在特异性脑电与大鼠觅药动机的形成及其行为的产生之间关联紧密。但由于海洛因诱导CPP模型大鼠在黑-白箱和白-黑穿梭过程中，不可避免地出现肌肉运动产生的肌电信号、电极摩擦噪声(其频带范围和实验大鼠主要节律脑电频带重叠)等干扰因素的影响[8]，常规脑电特征提取方法很难实现对实验大鼠个体觅药行为获得较高的实时分类准确率，从而无法对其强迫性觅药行为及时预警并对症干预。

本研究拟针对海洛因诱导CPP大鼠FrA区无线遥测脑电信号采取局部小波振幅阈值消噪预处理，再对消噪后脑电提取δ波、θ波、α波、β波，以此4种主要节律对应频率脑电小波系数的标准差值(线性特征)及样本熵值(非线性特征)作为融合特征向量，利用支持向量机(SVM)算法对成瘾大鼠个体不同行为状态的实时识别，结果得到较高的正确识别率，由此为海洛因诱导大鼠觅药行为的发生提供一个较为准确的检测方法。

1 材料与方法

1.1 海洛因诱导CPP大鼠模型制备及FrA区无线遥测场电位记录[4-6]

1.1.1 动物和frA区电极埋藏 清洁级 Wistar 大鼠( 南京市江宁区青龙山动物繁殖场，许可证号为 SCXK 苏 2001 -0001) 30只。在1%戊巴比妥钠(5 ml· kg-1)麻醉状态下行FrA区电极埋藏手术，左右侧FrA区定位：前卤+5.00 mm，矢状缝旁开 2.00 mm，深度2.5 mm，电极使用0.3 mm漆包镍铬线并用牙科自凝水泥固定，封闭切口后连续3 d注射青霉素以防手术感染。

1.1.2 海洛因诱导大鼠建模 随机选择并标记编码FrA电极埋藏的20只大鼠，作为海洛因诱导CPP组，放入 CPP 视频黑白箱内( JLBehv 条件性位置偏爱视频分析系统，上海吉量软件科技有限公司生产)，以白箱作为拌药箱并用活动挡板隔离黑、白箱，第1日上午9:00第1次在白箱内对其皮下注射海洛因0.5 mg/(kg·d)，第2日开始每天2次(上午9:00 1次，下午6:00 1次)，每日递增0.25 mg/kg,连续注射7 d，剩下10只大鼠(已施行FrA电极埋藏手术)作为手术对照组，同法注射等量生理盐水。海洛因诱导CPP大鼠注射7 d后，分别在停药后第1日和第2日上午做位置偏爱测定，让2组大鼠在黑-白箱自由行动，记录45 min黑和白箱行动轨迹和停留时间，做组间和海洛因诱导CPP组大鼠训练前后统计学分析。统计学Levene 法检验方差齐性，两独立样本比较用t检验，显著水准α=0.05，不满足检验条件时，采用秩检验。同一大鼠海洛因诱导前后拌药箱停留时间对比超过50%，可确立该大鼠个体海洛因诱导CPP制模成功。

1.1.3 2组大鼠FrA自发脑电测定 利用脑电生理无线遥测系统( BW - 200 型生理无线遥测系统，成都泰盟科技有限公司生产)及条件性位置偏爱视频系统每天2次(上午8:50—9:10，下午5:50—6:10)分别记录海洛因诱导大鼠自然戒断期(连续注射海洛因7 d后停止注射后24～48 h内)FrA区20 min脑电数据，并按01白-黑穿梭、02黑箱停留、03黑-白穿梭、04白箱停留标注其行为状态，每只大鼠总共记录4次共80 min数据。因为实验大鼠在4种行为状态之间多次切换且时间不等，为考虑每种行为状态数据长度的可比性，我们以2 s为脑电数据长度，截取经过后续1.2节中小波消噪后，每种行为状态对应脑电数据120段。最终本文分析的实验大鼠脑电数据形式为20×4×120，其中20代表戒断实验大鼠个数, 4代表其黑、白箱停留、白-黑穿梭、黑-白穿梭四种行为状态，120代表脑电样本组数。

1.2 基于特定肌电噪声区域的FrA区无线遥测脑电数据预处理方法

现阶段脑电消噪技术多针对于可放置噪声参考电极的多电极脑电原始信号且消噪效果较好，能用于单电极脑电消噪的技术主要有：EMD-ICA、TESWT自适应及小波变换方法[9]。小波变换因兼顾脑电信号的时域和频域信息在单电极脑电消噪中应用相对广泛[10,11]。加之，本课题组所研究的海洛因诱导大鼠自然戒断期脑电属于单电极脑电且肌电噪声相对较多。我们选用小波阈值消噪方法来预处理原始脑电信号。根据Mallat算法公式(1)可知，可对本实验采集的离散EEG单电极信号x(n)可进行有限层分解，

(1)

其中，L为分解层数，AL为低通逼近分量，Dj为不同尺度下的细节分量。因本实验脑电信号的采样率为500 Hz，为让小波不同尺度分解后系数能较好的对应脑电常见四节律δ波(0～3.9 Hz)、θ波(3.9～ 7.8 Hz)、α波(7.8～15.6 Hz)、β波(15.6～31.25 Hz)的频率范围及对应能量。我们选用db5小波对每段脑电数据进行6层分解时，则上述四个节律振幅依次分别可用A6、D6、D5、D4四个小波系数来表征。但这种消噪方法多是直接作用于整段脑电信号x(n)，导致面对本课题组所采集的其间包含和δ波、θ波、α波、β波主要节律频带范围接近而振幅远大于实验大鼠脑电的骨骼肌肌电的海洛因诱导大鼠自然戒断期脑电，消噪过程中因阈值的大小选取很难权衡而影响最终的消噪效果。针对上述肌电信号振幅远大于常规脑电振幅；所有消噪信号处理方法均可能对实验大鼠有效脑电信号有损耗等因素考虑。我们首先对每只戒断实验大鼠的80 m脑电数据按4种行为状态进行分类，然后计算每只大鼠每种行为状态所有脑电数据的平均振幅，最后得到20只实验戒断大鼠每种行为状态的平均振幅。若某只大鼠该种行为状态2 s的脑电数据平均振幅大于该平均值，则认为截取的该段数据含有肌电噪声信号；需对该段脑电信号首先采用db5小波进行6层分解，再对每一层分小波系数手动作硬阈值函数处理，最后再采用逆小波变换还原出消噪后的脑电信号；反之，则认为该2 s段原始脑电信号含肌电噪声干扰较少而无需再做消噪处理。因为从信号消噪处理原理的角度考虑，所有信号消噪处理方法均会在消除噪声的同时对原始脑电信号的有效成分也产生一定程度的衰减。最终将经过小波消噪预处理后采用逆小波还原出的某段2 s脑电信号和含肌电噪声干扰较少而未做消噪处理的2 s原始脑电信号汇总为某只大鼠该种行为状态下所有脑电数据。按1.1中方法预处理方法处理后得到实验大鼠每种行为状态下脑电样本数据120段，随机选取其中80段样本数据组成训练集数据(80×4=320段)，每种行为状态剩下的40段样本数据组成测试集数据(40×4=160段)。

1.3 海洛因CPP大鼠FrA区无线遥测脑电特征提取

海洛因CPP大鼠在动机形成及产生强迫性觅药行为时，其负责整合当前感觉信息、奖赏效价和机体需要的前额联络皮层区域脑电一定会呈现一定的特异性改变。信号标准差是信号幅值偏离平均值的均方根，文献[8]研究表明海洛因位置偏爱大鼠FrA区脑电中δ波、θ波、α波、β波四种主要节律和正常对照组相比，其能量值及含量百分比均出现不同程度的特异性改变，提示若能提取不同行为状态脑电信号的标准差以及对其进行小波分解后提取出与δ波、θ波、α波、β波四个主要节律频带对应的A6、D6、D5、D4四个小波系数的均方根作为四种不同行为状态下线性特征会对提高识别率有帮助。文献[7]研究表明海洛因CPP大鼠FrA区脑电样本熵值和正常对照组相比，也表现出较显著的改变，提示若能提取不同行为状态脑电信号的样本熵值作为四种不同行为状态下非线性特征也会对提高识别率有帮助。综上所述，经过这样融合处理后输入不同行为状态分类器的脑电特征量有4个线性和1个非线性特征共5个特征量。

1.4 支持向量机分类算法原理

支持向量机是一种能较好处理小样本监督分类，并具有很好的泛华能力的统计学习方法[11]。本研究选择的SVM分类模型为C-支持向量分类模型。其算法原理具体步骤如下：

(1)选择合适的核函数K(Xi，Xj)和惩罚参数C，构造并求解最优化问题

(2)

(3)

(4)

(3)构造决策函数：

(5)

公式(5)中K，ai*和b*分别代表支持向量个数，支持向量系数和分类阈值。本研究选用高斯径向基核函数作为核函数，其表达式为：

(6)

利用LibSVM工具箱，利用粒子群优化算法进行参数寻优，得到最佳惩罚因子c和核函数参数g。

2 结果

2.1 海洛因依赖(拌药箱位置偏爱)大鼠的确定

在大鼠45 min黑白箱穿梭时间中，选择前15 min穿梭轨迹做白箱(拌药箱)停留时间和停留时间百分比测定(后30 min大鼠有疲劳懈怠现象)，结果见表1。统计学显示，海洛因诱导CPP组大鼠停药1～2 d，其白箱停留时间及其百分比均较用药前明显增加(P<0.05，P<0.01)。与手术对照组(注射盐水)比较，上述指标更为显著(P<0.05，P<0.01)。上述结果表明本海洛因制模方法已成功使大鼠产生海洛因依赖。

Tab. 1 Analysis of staying time ( 15 min) of heroin-induced CPP rats in white box

2.2 基于特定肌电噪声区域的小波消噪预处理结果

为了便于解释消噪处理方法，我们特意挑选连续依次包含01白-黑穿梭、02黑箱停留、03黑-白穿梭、04白箱停留四个行为状态，总时长为60 s的4#海洛因位置偏爱诱导实验大鼠戒断期FrA区原始脑电(因实验过程中CPP实验大鼠在四种行为状态之间多次切换且时间长短严重不等)，做出脑电时域波形图如图1所示。从图1中(a)中可以看出实验大鼠四种行为状态中均有噪声污染带来的尖刺信号，尤以03黑-白穿梭状态最为明显(其振幅上下差值接近达到35 μV)。图1中(b)为对60 s原始脑电信号整体采用6层小波分解提取出对应δ波、θ波、α波、β波四个小波系数后对再采用硬阈值消噪方法处理，最后采用逆小波变换还原得出的全局阈值小波消噪后信号。与(a)中原始脑电信号相比，突兀的尖刺噪声信号基本消除，其中03黑-白穿梭状态的尖刺信号振幅下降明显(其振幅上下差值已削弱为接近20 μV左右)。图1中(c)为针对包含肌电噪声较多的原始脑电局部区域采用6层小波分解提取出对应δ波、θ波、α波、β波四个小波系数后对采用硬阈值消噪方法处理最后还原汇总得到的CPP组局部手工调整阈值小波消噪后脑电时域信号，与图1中(a)、(b)中脑电信号相比可看出，其实现了在尖刺噪声信号消除的同时原始脑电信号波形中的锐利有效成分得以保留。其中对03黑-白穿梭状态的尖刺信号最为明显，在基本保留了原始脑电信号中该段特征波形的同时将其尖刺信号振幅也已削弱为接近20 μV左右。

Fig. 1 Time domain waveform of CPP heroin-withdrawal rats

图2为与图1中三种情况下脑电时域波形图对应的频谱图，从图2(a)图可看出噪声频段主要位于高幅慢波的δ波(0～3.9 Hz)和θ波(3.9～7.8 Hz)频段范围内，其和振幅对应的强度值最高达到5 200左右。与原始脑电信号频谱(a)图相比，经过小波阈值消噪处理后的(b)、(c)频谱图均在横坐标为 0.5 HZ左右处有较大凹陷，最大振幅对应的强度值也从5 200左右削弱到4 200左右，说明小波阈值消噪方法对本研究所得到的原始脑电信号中所包含的肌电噪声消除非常有针对性。从图2(c)图可看出，针对特定肌电区域采用小波阈值消噪方法，除了在消除高幅慢波的δ波(0～3.9 Hz)和θ波(3.9～7.8 Hz)频段范围内的肌电噪声效果明显的同时保留了原始脑电中α波(7.8～15.6 Hz)、β波(15.6～31.25 Hz)频段有效信号的强度，从而为后续实验大鼠脑电特征提取及行为状态识别打下了良好的基础。

Fig. 2 Frequency domain waveform of CPP heroin-withdrawal rats

2.3 4#海洛因诱导CPP大鼠四种行为状态特征量的比较

盒形图是在1977年由美国的统计学家约翰.图基(John Tukey)发明的，能够有效帮助我们识别数据的特征[12]。它由上四分位数(Q3)、中位数和下四分位数(Q1)组成一个“带有隔间的盒子”。盒子的长度(IQR)等于上下四分位线之间的差值，最大观察值等于上四分位数和1.5倍的盒子长度之间的差值，最小观察值等于下四分位数和1.5倍的盒子长度之间的差值。一般上下四分位数和最大观察值和最小观察值之间有一条被称为“胡须”的延伸线。大于和小于最小观察值的数据点称为“离群点”，通常被单独标注以免因少数的“离群点”数据导致整体特性的偏移。对每个实验大鼠80 m的原始脑电数据做完针对特定区域肌电小波消噪处理后的脑电数据按行为状态分类各选取120个片段样本，每个片段长度为1 000个数据点(2 s)。提取每个样本脑电数据片段与δ波、θ波、α波、β波四个主要节律频带对应的A6、D6、D5、D4小波系数标准差(线性特征)及其样本熵值(非线性特征)，采用非参数检验来分析每只大鼠四种行为状态下δ波、θ波、α波、β波四个主要节律频带强度标准差值之间的差异。

图3为4#海洛因诱导CPP大鼠四种行为状态四种主要节律强度标准差值盒形图，显示4#实验大鼠额叶脑电数据在03黑-白穿梭行为状态下其δ波、θ波、α波、β波的强度标准差值和其它行为状态相比均存在显著差异，也即针对特定肌电噪声区域采用小波阈值消噪后得到的脑电数据的A6、D6、D5、D4小波系数标准差值(线性特征)适合作为海洛因诱导CPP大鼠强迫性觅药行为实时识别的特征值。但同时我们也看到四种主要节律在每种行为状态下均有少数的“离群点”，特别是图中(a)图可看出频带范围在0～3.9 HZ与δ波近似的“离群点”噪声数据相对较多，说明针对本研究的脑电消噪技术还有提高的空间。

Fig. 3 Standard difference box plot of main rhythm intensity from four behavioral states in 4# CPP heroin-withdrawal rats

图4为海洛因诱导CPP4#大鼠四种行为状态脑电样本熵值盒形图，显示4#实验大鼠额叶脑电数据在03黑-白穿梭行为状态下其样本熵值和其它行为状态相比均存在显著差异且明显小于其它三种行为状态。也即针对特定肌电噪声区域采用小波阈值消噪后得到的表征脑电数据序列复杂度的样本熵值(非线性特征)适合作为海洛因诱导CPP大鼠强迫性觅药行为实时识别的特征值。

Fig. 4 Sample entropy value box plot from four behavioral states in 4# CPP heroin-withdrawal rats

2.4 海洛因诱导CPP大鼠四种不同行为状态SVM分类准确率

表2为20只实验大鼠基于四个节律强度标准差及样本熵作为融合特征值到的01白-黑穿梭、02黑箱停留、03黑-白穿梭、04白箱停留四个行为状态的SVM分类准确率。选取单只实验大鼠四种行为状态额叶脑电各120样本片段(数据点为1 000，时长为2 s)，随机选取每种行为状态80个样本片段数据作为训练数据(总数据为80×4等于320个数据)，剩下的40个样本片段数据作为测试数据(总数据为40×4等于160个数据)，采用同样的算法计算出每只实验大鼠四种行为状态的分类准确率。从表中可以看出所有实验大鼠03黑-白穿梭状态的分类识别率较高，平均值达到83.88。

Tab. 2 Classification accuracy of four behavioral states in CPP heroin-withdrawal rats

3 讨论

海洛因药物加场景的条件性刺激，可诱导大鼠产生海洛因CPP，后者被广泛地运用于大鼠药物依赖模型的制作。本实验利用大鼠天然喜爱黑(暗)环境习性，将CPP视频白箱作为拌药箱，连续7 d递增剂量给予海洛因皮下注射并随即戒断，经视频观察大鼠黑-白箱穿梭轨迹和统计学分析，发现海洛因诱导CPP组大鼠注射药物7 d后戒断期，白箱停留时间和停留时间百分比显著增加，提示该组大鼠已产生了明显的白箱觅药行为和海洛因药物依赖。

海洛因诱导CPP大鼠突发的强迫性觅药行为与FrA区的脑电特征性改变之间关联紧密[13]。但诱导CPP大鼠在出现强迫性觅药行为时多伴随有与之相一致的肌电等干扰噪声[14]，而这些噪声的信号频带范围(0～8Hz)和实验大鼠强迫性觅药行为关联最紧密的特征脑电δ波(0～3.9 Hz)和θ波(3.9～7.8 Hz)频带大部重叠[5]，加之不同海洛因依赖实验大鼠个体脑电之间存在差异性等原因[15]，使得有关药物依赖大鼠强迫性觅药行为的脑电特征研究很难用于预警海洛因依赖患者个体的强迫性觅药行为[16]，及时对症干预，从而有效抑制并避免“复吸”。

本研究通过分析海洛因诱导CPP大鼠FrA区原始场电位信号，发现与穿梭觅药相一致的阵发性肌电噪声干扰信号其信号频带范围与高幅慢波的δ波(0～3.9 Hz)和θ波(3.9～7.8 Hz)频段范围接近，但振幅远大于实验大鼠与强迫性觅药行为关联紧密的脑电有效信号振幅。针对肌电噪声干扰信号的上述特点，我们对包含特定噪声较多而表现出振幅过大区域针对性采用小波阈值消噪方法对原始脑电信号进行消噪处理。从海洛因诱导CPP戒断大鼠时域波形图1和频域波形图2中可以看出，本文所采用的消噪方法在取得很好的消噪效果的同时最大程度的保留了原始信号中的特征脑电有效成分。海洛因诱导CPP大鼠个体穿梭觅药行为相关脑电消噪和特征值的实时提取，有助于海洛因诱导CPP大鼠真实脑电信号的识别，提高脑电信号与成瘾大鼠个体觅药行为之间的相关性。

脑电信号处理的方法从原理上可分为线性分析和非线性分析二大类[17]。有研究表明，与海洛因诱导CPP大鼠脑电中和δ波、θ波、α波、β波四个主要节律对应的小波系数标准差所表征的线性特征变化，以及表征脑电序列信号复杂性的样本熵值标准差所代表的非线性变化均与实验大鼠的觅药行为关联紧密[13，14]。所以，本研究中我们对经过前述针对特定区域采用小波阈值消噪得到的实验大鼠脑电，首先按行为状态进行分类各自选取时长为2 s的120个脑电样本片段，再提取其对应四种主要节律强度标准差值和样本熵值。从实验大鼠4种行为状态和4种主要节律强度标准差值盒形图3，脑电样本熵值盒形图4可看出，上述4个线性特征和1个非线性特征均非常适合作为海洛因诱导CPP大鼠强迫性觅药行为实时识别的特征值，融合后的特征值因既考虑了(FrA)原始脑电的线性特征又兼顾了其非线性特征，从而具有更好的分类识别率。从表1海洛因诱导CPP大鼠4种行为状态分类准确率可知：20只实验大鼠个体的四种行为状态分类准确率均值都达到在80%上下浮动，其中4#实验大鼠的03黑-白穿梭达到了83.88%的分类识别率，基本实现了海洛因CPP大鼠行为状态的自动分类与预警。

综上所述，针对海洛因诱导CPP大鼠个体FrA脑电原始数据，采用小波阈值消噪预处理去除特定肌电噪声，再选取所有实验大鼠4种行为状态FrA脑电各120样本片段提取其δ波、θ波、α波、β波，同时对这些节律小波系数标准差值和样本熵值组成融合特征向量，输入支持向量机，每只实验大鼠个体的4种行为状态分类准确率均值均达到预想的结果，基本实现了对海洛因成瘾大鼠个体强迫性觅药行为的实时识别，不失为一种检测海洛因大鼠觅药行为发动、发生的有效方法。籍此思路，后续开展更大规模样本数的探索研究或者直接用无线脑电监测设备获得海洛因依赖患者脑电数据并试图探索其与觅药紧密相关的行为状态之间的关系并有效识别,对于应用于海洛因依赖患者的临床观察和预防觅药行为的产生，具有一定的借鉴和研究价值。