目前，帕金森病发生率在逐渐上升，严重影响着患者生活质量，社会诊疗负担不断加重。然而，该病的早期监测手段有限，很难及时干预。为了发现其标志物，本文对服药前后帕金森病患者及健康人的32通道静息态脑电数据进行分频段研究。首先利用动态加权符号互信息计算各通道脑电信号间相关性，再通过k均值聚类实现信号关联矩阵的分类，最后得到脑电信号关联状态。通过统计分析发现，在Beta频段（P = 0.034）与Gamma频段（P = 0.010）各有一个脑电信号关联状态可显著区分未服药帕金森病患者与健康人。这表明未服药帕金森病患者与健康人的静息态脑电各通道信号相关性差异有统计学意义。而在服药与未服药帕金森病患者、服药帕金森病患者与健康人之间，其关联状态差异均没有统计学意义。这可为帕金森病的临床诊断提供一种参考。
引用本文： 丁昊, 吴进辉, 唐旭东, 余江南, 陈轩恒, 吴占雄. 基于动态加权符号互信息与k均值聚类的帕金森病患者静息脑电关联状态识别. 生物医学工程学杂志, 2023, 40(1): 20-26. doi: 10.7507/1001-5515.202211002 复制
0 引言
帕金森病（Parkinson’s disease，PD）是一种常见的慢性行进性神经系统疾病，多见于老年人，给患者及其家庭带来极大的痛苦与严重的经济负担。它主要是由黑质多巴胺能神经元的异常变化引起的[1-2]。帕金森病主要表现为运动缓慢、姿势不稳、静止性震颤，以及认知、睡眠、精神与自主神经功能障碍等非运动症状[3-5]。帕金森病的治疗手段主要有药物、手术、神经康复等，也可以通过增加多巴胺和神经元保护以延缓疾病进展。然而，其确切发病机制目前尚不清楚，也没有药物可以逆转病程以实现根治[6]。发现简单易用的帕金森病生物学标志，有助于及早发现、干预疾病，将为患者带来福音。
脑电信号（electroencephalogram，EEG）利用非侵入式电极记录脑神经元活动情况，能够反映大脑皮层的活动规律，可为帕金森病诊断提供生物学标志[7-8]。目前已经有大量研究表明，帕金森病患者相比于健康人在脑电图上存在异常。多巴胺神经元的变性衰亡会导致基底节-丘脑-皮层回路的同步模式明显异常，进一步影响帕金森病患者的脑网络信息编码[9]。电生理学研究表明，帕金森病表现为基底节-皮层环路、皮层-皮层环路的病理振荡，其中基底节同步活动与皮层的活动呈线性相关[10]。通过对脑电信号进行分析，可以发现基底节-丘脑-皮层回路活动改变的标记，为了解其大脑皮层异常动态变化提供新的途径[9]。帕金森病患者的枕部Alpha与Theta频段脑电活动信号减弱，视觉空间障碍会导致枕部右侧脑电信号异常[11]。未服药帕金森病患者的Beta频段脑电信号波形在感觉运动电极上的振荡比服药后有更大的锐度不对称和陡度不对称性[12]。处于静息状态时，脑电图显示帕金森病患者与健康人对照组相比在Theta和Alpha频段上功率增加，Beta和Gamma频段上功率降低[13]。在脑电信号频谱图中，帕金森病患者的Beta频段波振幅比健康人小[14]，而Alpha波段脑电信号振幅比健康人大，且其频率减慢[15]。Zhang等[16]通过将时频分析与深度学习相结合，应用可调Q因子小波变换深度残余收缩网络来对脑电数据进行分类，帕金森病预测的准确率达到99.92%。在对高分辨脑电数据源定位后，Gerard等[17]采用动态锁相方法来评估帕金森病患者功能连接性，发现患者额区与枕颞区显示出异常强的Theta频带连接性。基于脑电数据5个频带的相干性，Conti等[18]通过计算连接矩阵及分析网络拓扑参数，发现帕金森病患者在Alpha、Beta频带连通性较低，而在Gamma频带连接性增强。此外，各通道脑电信号并非独立不相关，它们之间相互影响，可能存在着某种非线性依赖关系，蕴含着可以作为检验帕金森病的生物学标志。为此，加权符号互信息算法（weighted symbolic mutual information，wSMI）可以用于探索脑动态功能连接性及病理条件下的脑网络潜在变化（如阿尔茨海默病、多发性硬化症、精神分裂症和社交焦虑障碍）[19]。有研究使用wSMI与聚类来研究脑连接的动态过程，并对比分析了精神分裂症患者与健康人的丘脑区域网络效率[20]。也有人采用wSMI研究脑电信号Delta、Theta频带的脑皮质功能连接性[21]。
基于以上思考，本文假设帕金森病患者与正常健康人的静息态脑电信号在各个频段内可能存在关联状态差异。在计算静息态脑电信号关联状态后，通过比较帕金森病患者与健康人的各通道脑电信号间的相关性来发现生物学标志。脑电信号预处理后，利用wSMI计算通道间相关性以得到动态连接矩阵，再利用k均值聚类进行关联状态识别，并计算出不同对象的关联状态发生比例，最后统计发现与帕金森病显著相关的功能连接状态。这里的状态是指不同通道脑电信号之间存在的神经活动同步性。
1 方法
1.1 脑电数据
本文所采用的脑电数据为加州大学圣地亚哥分校帕金森病患者静息状态脑电开放数据[22]。数据集包括15名帕金森病患者［8名女性，7名男性，年龄为（63.27±8.20岁）］和15名正常健康人［9名女性，6名男性，年龄为（63.50±9.67岁）］，所有人均为右利手。在采集脑电数据时，患者继续药物治疗方案，治疗药物主要为能影响多巴胺能的左旋多巴等。通过平衡患者用药和停药的时间，分别收集患者按时用药的脑电信号（PD-on）以及患者至少停止用药12 h（PD-off）的脑电信号。同时，收集健康人（healthy control，HC）的脑电信号作为对比分析。脑电信号采集对象信息如表1所示。脑电信号采用32通道的BioSemi ActiveTwo系统采集，并有额外的八个电极放置在眼睛的外侧和下方，以监测眼睛的眨眼和运动。采集时所有人均处于舒适的静息状态，采样频率为512 Hz，静息数据连续记录至少3 nub。这里值得注意的是，PD-off与PD-on脑电信号采集于同一组PD患者。
表1
研究对象信息
Table1.
Demographic information of the subjects
1.2 脑电信号预处理
将脑电信号导入到EEGLAB MATLAB工具箱中进行预处理，具体步骤如下：① 坐标定位：对电极进行坐标定位，只对前32个通道采集的信号进行分析，去除后8个不反映脑电信号的电极（EXG1，EXG2，…，EXG8）。② 重定义参考：取每个电极的平均值，应用全部脑电数据的共同平均值进行重定义参考。③ 滤波：使用四阶巴特沃斯数字滤波器带通滤波，分别得到每个对象脑电信号的Delta（0.5～4 Hz）、Theta（4～8 Hz）、Alpha（8～13 Hz）、Beta（13～30 Hz）、Gamma（30～50 Hz）及高频High-frequency（50～100 Hz）等六个频带区间，过滤掉其他过高频噪声信号和过低频的漂移伪迹。④ 独立成分分析：通过独立成分分析（independent component analysis，ICA），去除眼电等生理伪影的干扰，得到纯净的脑电信号。⑤ 截段提取：考虑到每个人在仪器测量的开始和结束可能会有心理状态的波动，未能立刻进入静息状态，同时各个对象所测得的脑电信号时间及数据序列长度有所不同，为了得到较为稳定的脑电信息，以脑电信号时间序列中心为基准向前后两侧各取80 s，剔除了数据最前段和最后段，只保留中间有效部分，最终得到160 s的数据，使所有对象的脑电信号时间序列长度均相等。
1.3 动态wSMI分析
本文采用wSMI来评估脑电信号随机非线性联合波动的程度，度量时间序列间的共享信息[21]。同时，该算法能够发现信号增加或减少的定性“符号”模式，在不需要对交互类型进行假设的情况下能对信号的熵进行快速和稳健估计。加权符号互信息通过计算不同时间序列间同时出现的每一对符号的联合概率来估计两个信号之间共享符号的互信息，并通过忽略相同或相反符号信号的共同出现以解决容积传导的问题[19]，其权值计算能够避免来自共同来源的脑电信号间的伪相关性，使得到的结论更为严格。这是一种检测非线性耦合的有效方法[23]。其计算公式如下：
 |
其中
、
为时间序列，
为
中元素，
为
中元素。
为构成符号的元素个数，
为加权函数，
为概率函数，
为自然对数函数。
首先根据预先定义的连续采样数据的振幅趋势将脑电信号转变为一系列离散符号。通过估计每对符号的联合概率来计算脑电信号序列间的加权符号互信息值。对于一对相同的符号，权值设置为零，表示这对符号可以由唯一的共同源引出[9]。为了分析脑电信号时间序列的动态特性，不是一次性对整个时间序列进行加权符号互信息计算，而是通过施加时间窗口，每次取一段脑电信号进行处理计算。随着时间窗口的滑动，可以得到一系列加权符号互信息计算结果，构成三维矩阵。这样可以更好地体现出脑状态的时间波动性。本文里，wSMI采用Python语言实现，底层调用MNE库（MEG、EEG analysis & visualization，用来处理EEG/MEG等信号的开源程序，
1.4 k均值聚类
k均值聚类，是一种经典的高效快速聚类方法，其计算过程是把多个样本分为k个簇，按照距离中心点最“小”的原则，使得各簇内的样本点具有较高的相似性，而划为不同簇的样本相似性较低。相似度是依据一个簇中样本的平均值来进行计算的。每次迭代后重新计算聚类中心，重新划分类别，由此不断迭代重复，直至聚类中心值趋于稳定[26]。
将加权符号互信息得到的关联矩阵集合按照每个频段分别聚类。这里聚类数目取k = 6[20]。对每个矩阵样本计算它到聚类中心的距离并将它分到距离最小的类中，然后针对每个矩阵样本类别重新计算聚类中心。兼顾速度和效率，这里设置重复迭代上限为200次。最后得出所有频段的聚类中心、样本所属类别、概率分布和转移概率等数据。最后得到的关联矩阵如图1所示。
图1
脑电信号关联矩阵
Figure1.
Correlation matrix between each pair of EEG channels
1.5 统计分析
不考虑其他因素的影响，把脑电数据分为健康人HC、未服药帕金森病患者PD-off、服药帕金森病患者PD-on三类（每组15人，共三组），作为ANOVA统计分组变量。在验证所得的样本数据符合正态分布后，计算HC、PD-off、PD-on各组对象在不同频段关联状态占比，之后对这些样本进行ANOVA统计分析。在发现不同组间差异有统计学意义后，采用最小显著差别方法（least significant difference，LSD）进行校验。检验水准为0.05。本论文所实现的算法源代码与分析结果数据可从
2 结果
图2展示了不同频段静息态脑电信号关联矩阵聚类结果，每个频带共聚类为6个关联状态。本论文采用状态所占比例作为指标来区分不同分组，探究将静息态脑电信号关联状态作为帕金森病生物学标志的可能性。健康人、服药帕金森病患者、未服药帕金森病患者三个分组的脑电信号关联状态占比数据如表2所示。
图2
对三组对象采用加权符号互信息及k均值聚类获取静息态脑电信号关联状态
每个频段聚类得到6个关联状态，采用关联状态占比作为指标来区分不同组
Figure2. Correlation matrices extracted from EEG signals using wSMI and k-means clustering in three groupssix correlation estimates were obtained for each frequency band. The correlation status proportion was used as an index to distinguish different groups
表2
三组对象静息态脑电信号关联状态占比（
±s）
Table2.
Correlation status proportion of HC, PD-off, and PD-on groups (
±s)
ANOVA统计分析之后，发现HC与PD-off组间Beta频带State1占比差异有统计学意义（P = 0.034，见表3），同时这两组间Gamma频带State5占比差异也有统计学意义（P = 0.010，见表4）。这说明在帕金森病患者与健康人之间，静息态脑电信号关联性在Beta和Gamma频带显著不同。
表3
Beta频带关联状态1占比的显著性水平（LSD校正）
Table3.
Significance level of the proportion of state1 in Beta band (LSD test)
表4
Gamma频带关联状态State5占比的显著性水平（LSD校正）
Table4.
Significance level of the proportion of state5 in Gamma band (LSD test)
3 讨论与结论
本论文利用动态wSMI及k均值聚类实现了帕金森病患者静息态脑电信号关联性分析。发现在Beta段范围内，健康人与未服药帕金森病患者状态State1占比差异有统计学意义（P = 0.034），而未服药与服药帕金森病患者之间State1占比差异没有统计学意义（P = 0.437），健康人与服药患者间State1占比差异也没有统计学意义（P = 0.169）。在Gamma频带，健康人与未服药帕金森病患者间State5占比差异有统计学意义（P = 0.010），而未服药与服药患者间State5占比差异没有统计学意义（P = 0.560），健康人与服药患者间State5占比差异也没有统计学意义（P=0.053）。
考虑到k均值聚类可能存在误差，本文进行了多次重复聚类计算，均得到同样的结果。同时，观察到Delta频带的状态图与其他频段有明显的不同（见图2），猜想在Delta频带通道相关性可能具有与其他频带完全不同的独特性质。此外，还发现帕金森病患者的脑电波形在后半段有更大的波动和陡峭性，而健康人的则平坦稳定，这与之前的研究结果是一致的[12]。在k均值聚类时，Delta与High-frequency频段具有比较大的波动性，有可能是频段过窄或过宽的原因。
另外，本文也验证了当聚类数目分别设为4、5、7时，并没有发现差异具有统计学意义的脑网络关联状态，而只有在聚类数为6时发现两个有显著差异的状态。这也验证了Fu等[20]的工作。
同时，分析本文研究方法可能产生误差的原因：① 独立成分分析去除生理伪影。采用独立成分分析去除生理伪影时，对生理伪影的判断可能会有一定的误差。② 高频带的噪声干扰会对互信息计算有所影响。③ 采集静息态数据时对象可能存在情绪或认知活动波动，例如患者脑神经功能的缺陷使得对象不能长时间集中精力等。
在Beta频带State1和Gamma频带State5占比上，尽管未服药帕金森病患者与服药帕金森病患者之间、服药帕金森病患者与健康人之间差异没有统计学意义，但还是可以观察到患者服药后其状态占比的统计均值介于未服药和健康人之间，说明服药能够缩小这两者状态的差距，因此服药对于病变的程度还是有所改善的。关于帕金森病患者服药与未服药的差异，还可以在此基础上进一步深入研究。同时，在本研究采用32通道脑电数据得出了上述结论，如果基于更多采集通道（如64、128等）的脑电信号数据，脑电信号关联状态的了解与分析可以更为透彻与精确。
wSMI能够评估脑电信号间非随机联合波动的程度，度量其共享信息值。相较于其他相关性算法，加权符号互信息主要有三个优点：① 能够寻找信号增加或减少的定性或“符号”模式，从而快速和准确地估计信号的熵，符号的转换取决于符号的长度和时间间隔；② 提供了一种检测非线性耦合的有效方法，并且不需要对相互作用的类型进行假设，提高了结果的真实性；③ 权重计算舍弃了脑电信号之间的伪相关性[21, 23]。这些特点使得wSMI非常合适脑电信号关联状态的识别。
总之，本研究采取动态加权符号互信息方法计算脑电信号间关联性，通过聚类评估帕金森病患者脑电信号关联状态的改变模式，从而发现帕金森病的生物学标志。这为帕金森病患者的病情诊断与筛查提供了一种参考方法。
重要声明
利益冲突声明：本文全体作者均声明不存在利益冲突。
作者贡献声明：丁昊负责数据分析和论文写作，吴进辉与唐旭东负责脑电数据预处理，余江南与陈轩恒负责算法编程，吴占雄负责试验设计与文字修改。
0 引言
帕金森病（Parkinson’s disease，PD）是一种常见的慢性行进性神经系统疾病，多见于老年人，给患者及其家庭带来极大的痛苦与严重的经济负担。它主要是由黑质多巴胺能神经元的异常变化引起的[1-2]。帕金森病主要表现为运动缓慢、姿势不稳、静止性震颤，以及认知、睡眠、精神与自主神经功能障碍等非运动症状[3-5]。帕金森病的治疗手段主要有药物、手术、神经康复等，也可以通过增加多巴胺和神经元保护以延缓疾病进展。然而，其确切发病机制目前尚不清楚，也没有药物可以逆转病程以实现根治[6]。发现简单易用的帕金森病生物学标志，有助于及早发现、干预疾病，将为患者带来福音。
脑电信号（electroencephalogram，EEG）利用非侵入式电极记录脑神经元活动情况，能够反映大脑皮层的活动规律，可为帕金森病诊断提供生物学标志[7-8]。目前已经有大量研究表明，帕金森病患者相比于健康人在脑电图上存在异常。多巴胺神经元的变性衰亡会导致基底节-丘脑-皮层回路的同步模式明显异常，进一步影响帕金森病患者的脑网络信息编码[9]。电生理学研究表明，帕金森病表现为基底节-皮层环路、皮层-皮层环路的病理振荡，其中基底节同步活动与皮层的活动呈线性相关[10]。通过对脑电信号进行分析，可以发现基底节-丘脑-皮层回路活动改变的标记，为了解其大脑皮层异常动态变化提供新的途径[9]。帕金森病患者的枕部Alpha与Theta频段脑电活动信号减弱，视觉空间障碍会导致枕部右侧脑电信号异常[11]。未服药帕金森病患者的Beta频段脑电信号波形在感觉运动电极上的振荡比服药后有更大的锐度不对称和陡度不对称性[12]。处于静息状态时，脑电图显示帕金森病患者与健康人对照组相比在Theta和Alpha频段上功率增加，Beta和Gamma频段上功率降低[13]。在脑电信号频谱图中，帕金森病患者的Beta频段波振幅比健康人小[14]，而Alpha波段脑电信号振幅比健康人大，且其频率减慢[15]。Zhang等[16]通过将时频分析与深度学习相结合，应用可调Q因子小波变换深度残余收缩网络来对脑电数据进行分类，帕金森病预测的准确率达到99.92%。在对高分辨脑电数据源定位后，Gerard等[17]采用动态锁相方法来评估帕金森病患者功能连接性，发现患者额区与枕颞区显示出异常强的Theta频带连接性。基于脑电数据5个频带的相干性，Conti等[18]通过计算连接矩阵及分析网络拓扑参数，发现帕金森病患者在Alpha、Beta频带连通性较低，而在Gamma频带连接性增强。此外，各通道脑电信号并非独立不相关，它们之间相互影响，可能存在着某种非线性依赖关系，蕴含着可以作为检验帕金森病的生物学标志。为此，加权符号互信息算法（weighted symbolic mutual information，wSMI）可以用于探索脑动态功能连接性及病理条件下的脑网络潜在变化（如阿尔茨海默病、多发性硬化症、精神分裂症和社交焦虑障碍）[19]。有研究使用wSMI与聚类来研究脑连接的动态过程，并对比分析了精神分裂症患者与健康人的丘脑区域网络效率[20]。也有人采用wSMI研究脑电信号Delta、Theta频带的脑皮质功能连接性[21]。
基于以上思考，本文假设帕金森病患者与正常健康人的静息态脑电信号在各个频段内可能存在关联状态差异。在计算静息态脑电信号关联状态后，通过比较帕金森病患者与健康人的各通道脑电信号间的相关性来发现生物学标志。脑电信号预处理后，利用wSMI计算通道间相关性以得到动态连接矩阵，再利用k均值聚类进行关联状态识别，并计算出不同对象的关联状态发生比例，最后统计发现与帕金森病显著相关的功能连接状态。这里的状态是指不同通道脑电信号之间存在的神经活动同步性。
1 方法
1.1 脑电数据
本文所采用的脑电数据为加州大学圣地亚哥分校帕金森病患者静息状态脑电开放数据[22]。数据集包括15名帕金森病患者［8名女性，7名男性，年龄为（63.27±8.20岁）］和15名正常健康人［9名女性，6名男性，年龄为（63.50±9.67岁）］，所有人均为右利手。在采集脑电数据时，患者继续药物治疗方案，治疗药物主要为能影响多巴胺能的左旋多巴等。通过平衡患者用药和停药的时间，分别收集患者按时用药的脑电信号（PD-on）以及患者至少停止用药12 h（PD-off）的脑电信号。同时，收集健康人（healthy control，HC）的脑电信号作为对比分析。脑电信号采集对象信息如表1所示。脑电信号采用32通道的BioSemi ActiveTwo系统采集，并有额外的八个电极放置在眼睛的外侧和下方，以监测眼睛的眨眼和运动。采集时所有人均处于舒适的静息状态，采样频率为512 Hz，静息数据连续记录至少3 nub。这里值得注意的是，PD-off与PD-on脑电信号采集于同一组PD患者。
表1
研究对象信息
Table1.
Demographic information of the subjects
1.2 脑电信号预处理
将脑电信号导入到EEGLAB MATLAB工具箱中进行预处理，具体步骤如下：① 坐标定位：对电极进行坐标定位，只对前32个通道采集的信号进行分析，去除后8个不反映脑电信号的电极（EXG1，EXG2，…，EXG8）。② 重定义参考：取每个电极的平均值，应用全部脑电数据的共同平均值进行重定义参考。③ 滤波：使用四阶巴特沃斯数字滤波器带通滤波，分别得到每个对象脑电信号的Delta（0.5～4 Hz）、Theta（4～8 Hz）、Alpha（8～13 Hz）、Beta（13～30 Hz）、Gamma（30～50 Hz）及高频High-frequency（50～100 Hz）等六个频带区间，过滤掉其他过高频噪声信号和过低频的漂移伪迹。④ 独立成分分析：通过独立成分分析（independent component analysis，ICA），去除眼电等生理伪影的干扰，得到纯净的脑电信号。⑤ 截段提取：考虑到每个人在仪器测量的开始和结束可能会有心理状态的波动，未能立刻进入静息状态，同时各个对象所测得的脑电信号时间及数据序列长度有所不同，为了得到较为稳定的脑电信息，以脑电信号时间序列中心为基准向前后两侧各取80 s，剔除了数据最前段和最后段，只保留中间有效部分，最终得到160 s的数据，使所有对象的脑电信号时间序列长度均相等。
1.3 动态wSMI分析
本文采用wSMI来评估脑电信号随机非线性联合波动的程度，度量时间序列间的共享信息[21]。同时，该算法能够发现信号增加或减少的定性“符号”模式，在不需要对交互类型进行假设的情况下能对信号的熵进行快速和稳健估计。加权符号互信息通过计算不同时间序列间同时出现的每一对符号的联合概率来估计两个信号之间共享符号的互信息，并通过忽略相同或相反符号信号的共同出现以解决容积传导的问题[19]，其权值计算能够避免来自共同来源的脑电信号间的伪相关性，使得到的结论更为严格。这是一种检测非线性耦合的有效方法[23]。其计算公式如下：
|
其中、为时间序列，为中元素，为中元素。为构成符号的元素个数，为加权函数，为概率函数，为自然对数函数。
首先根据预先定义的连续采样数据的振幅趋势将脑电信号转变为一系列离散符号。通过估计每对符号的联合概率来计算脑电信号序列间的加权符号互信息值。对于一对相同的符号，权值设置为零，表示这对符号可以由唯一的共同源引出[9]。为了分析脑电信号时间序列的动态特性，不是一次性对整个时间序列进行加权符号互信息计算，而是通过施加时间窗口，每次取一段脑电信号进行处理计算。随着时间窗口的滑动，可以得到一系列加权符号互信息计算结果，构成三维矩阵。这样可以更好地体现出脑状态的时间波动性。本文里，wSMI采用Python语言实现，底层调用MNE库（MEG、EEG analysis & visualization，用来处理EEG/MEG等信号的开源程序，
1.4 k均值聚类
k均值聚类，是一种经典的高效快速聚类方法，其计算过程是把多个样本分为k个簇，按照距离中心点最“小”的原则，使得各簇内的样本点具有较高的相似性，而划为不同簇的样本相似性较低。相似度是依据一个簇中样本的平均值来进行计算的。每次迭代后重新计算聚类中心，重新划分类别，由此不断迭代重复，直至聚类中心值趋于稳定[26]。
将加权符号互信息得到的关联矩阵集合按照每个频段分别聚类。这里聚类数目取k = 6[20]。对每个矩阵样本计算它到聚类中心的距离并将它分到距离最小的类中，然后针对每个矩阵样本类别重新计算聚类中心。兼顾速度和效率，这里设置重复迭代上限为200次。最后得出所有频段的聚类中心、样本所属类别、概率分布和转移概率等数据。最后得到的关联矩阵如图1所示。
图1
脑电信号关联矩阵
Figure1.
Correlation matrix between each pair of EEG channels
1.5 统计分析
不考虑其他因素的影响，把脑电数据分为健康人HC、未服药帕金森病患者PD-off、服药帕金森病患者PD-on三类（每组15人，共三组），作为ANOVA统计分组变量。在验证所得的样本数据符合正态分布后，计算HC、PD-off、PD-on各组对象在不同频段关联状态占比，之后对这些样本进行ANOVA统计分析。在发现不同组间差异有统计学意义后，采用最小显著差别方法（least significant difference，LSD）进行校验。检验水准为0.05。本论文所实现的算法源代码与分析结果数据可从
2 结果
图2展示了不同频段静息态脑电信号关联矩阵聚类结果，每个频带共聚类为6个关联状态。本论文采用状态所占比例作为指标来区分不同分组，探究将静息态脑电信号关联状态作为帕金森病生物学标志的可能性。健康人、服药帕金森病患者、未服药帕金森病患者三个分组的脑电信号关联状态占比数据如表2所示。
图2
对三组对象采用加权符号互信息及k均值聚类获取静息态脑电信号关联状态
每个频段聚类得到6个关联状态，采用关联状态占比作为指标来区分不同组
Figure2. Correlation matrices extracted from EEG signals using wSMI and k-means clustering in three groupssix correlation estimates were obtained for each frequency band. The correlation status proportion was used as an index to distinguish different groups
表2
三组对象静息态脑电信号关联状态占比（±s）
Table2.
Correlation status proportion of HC, PD-off, and PD-on groups (±s)
ANOVA统计分析之后，发现HC与PD-off组间Beta频带State1占比差异有统计学意义（P = 0.034，见表3），同时这两组间Gamma频带State5占比差异也有统计学意义（P = 0.010，见表4）。这说明在帕金森病患者与健康人之间，静息态脑电信号关联性在Beta和Gamma频带显著不同。
表3
Beta频带关联状态1占比的显著性水平（LSD校正）
Table3.
Significance level of the proportion of state1 in Beta band (LSD test)
表4
Gamma频带关联状态State5占比的显著性水平（LSD校正）
Table4.
Significance level of the proportion of state5 in Gamma band (LSD test)
3 讨论与结论
本论文利用动态wSMI及k均值聚类实现了帕金森病患者静息态脑电信号关联性分析。发现在Beta段范围内，健康人与未服药帕金森病患者状态State1占比差异有统计学意义（P = 0.034），而未服药与服药帕金森病患者之间State1占比差异没有统计学意义（P = 0.437），健康人与服药患者间State1占比差异也没有统计学意义（P = 0.169）。在Gamma频带，健康人与未服药帕金森病患者间State5占比差异有统计学意义（P = 0.010），而未服药与服药患者间State5占比差异没有统计学意义（P = 0.560），健康人与服药患者间State5占比差异也没有统计学意义（P=0.053）。
考虑到k均值聚类可能存在误差，本文进行了多次重复聚类计算，均得到同样的结果。同时，观察到Delta频带的状态图与其他频段有明显的不同（见图2），猜想在Delta频带通道相关性可能具有与其他频带完全不同的独特性质。此外，还发现帕金森病患者的脑电波形在后半段有更大的波动和陡峭性，而健康人的则平坦稳定，这与之前的研究结果是一致的[12]。在k均值聚类时，Delta与High-frequency频段具有比较大的波动性，有可能是频段过窄或过宽的原因。
另外，本文也验证了当聚类数目分别设为4、5、7时，并没有发现差异具有统计学意义的脑网络关联状态，而只有在聚类数为6时发现两个有显著差异的状态。这也验证了Fu等[20]的工作。
同时，分析本文研究方法可能产生误差的原因：① 独立成分分析去除生理伪影。采用独立成分分析去除生理伪影时，对生理伪影的判断可能会有一定的误差。② 高频带的噪声干扰会对互信息计算有所影响。③ 采集静息态数据时对象可能存在情绪或认知活动波动，例如患者脑神经功能的缺陷使得对象不能长时间集中精力等。
在Beta频带State1和Gamma频带State5占比上，尽管未服药帕金森病患者与服药帕金森病患者之间、服药帕金森病患者与健康人之间差异没有统计学意义，但还是可以观察到患者服药后其状态占比的统计均值介于未服药和健康人之间，说明服药能够缩小这两者状态的差距，因此服药对于病变的程度还是有所改善的。关于帕金森病患者服药与未服药的差异，还可以在此基础上进一步深入研究。同时，在本研究采用32通道脑电数据得出了上述结论，如果基于更多采集通道（如64、128等）的脑电信号数据，脑电信号关联状态的了解与分析可以更为透彻与精确。
wSMI能够评估脑电信号间非随机联合波动的程度，度量其共享信息值。相较于其他相关性算法，加权符号互信息主要有三个优点：① 能够寻找信号增加或减少的定性或“符号”模式，从而快速和准确地估计信号的熵，符号的转换取决于符号的长度和时间间隔；② 提供了一种检测非线性耦合的有效方法，并且不需要对相互作用的类型进行假设，提高了结果的真实性；③ 权重计算舍弃了脑电信号之间的伪相关性[21, 23]。这些特点使得wSMI非常合适脑电信号关联状态的识别。
总之，本研究采取动态加权符号互信息方法计算脑电信号间关联性，通过聚类评估帕金森病患者脑电信号关联状态的改变模式，从而发现帕金森病的生物学标志。这为帕金森病患者的病情诊断与筛查提供了一种参考方法。
重要声明
利益冲突声明：本文全体作者均声明不存在利益冲突。
作者贡献声明：丁昊负责数据分析和论文写作，吴进辉与唐旭东负责脑电数据预处理，余江南与陈轩恒负责算法编程，吴占雄负责试验设计与文字修改。
