基于局部差异最小生成树脑网络的抑郁症分类研究

2020-07-22李舒婷王昕璨赵云芃陈俊杰

太原理工大学学报 2020年4期

李舒婷，李瑶，王昕璨，赵云芃，陈俊杰

(太原理工大学 a.信息与计算机学院，b.艺术学院，太原 030024)

复杂脑网络的分析与研究是近几年来神经精神疾病领域的研究热点。作为复杂网络理论在神经认知科学的具体应用，复杂脑网络在了解有关神经精神疾病的发病机理方面起到了很重要的作用，也为相关影像学标志物以及脑疾病临床的诊断和评价贡献了潜在价值和新的方法。

尽管目前复杂脑网络领域的研究受到众多学者广泛关注，但是由于网络分析技术的不成熟、特征选择及分类模型构建方法的多样性等，这一领域仍然存在着诸多亟待解决的问题。比如，传统复杂网络分析方法在不同组间进行对比的标准化问题[1]。标准化指的是网络的规模，包括网络大小(即节点的数目)、网络稀疏度(即存在连接的占比)以及平均度(即每个节点的连接数)等。需要注意的是，目前该领域内对复杂脑网络的稀疏度尚未形成统一的结论。基于此，有研究人员提出利用最小生成树(minimum spanning tree,MST)方法来进行全脑范围的复杂脑网络的构建[2-3]。这一无偏方法可以在保证所有节点均连入网络的前提下，网络中的边被最大化的精简，网络结构得到极大的简化，而网络核心框架仍可以被保留。因此，最小生成树网络在保证网络连接的同时，尽可能地保持较高的连接强度。通过这一方法，在一定程度上解决了网络稀疏度等因素配置对网络结构的影响，已被广泛应用于脑疾病的研究中。

然而，在全脑构建MST的方法仍存在一些局限。在这一方法中，不考虑所有的连接以防止聚类，从而获得固定的网络大小和密度，以使网络可以比较。因此，MST可能低估了其他有趣信息的重要性，如低权重连接、脑网络中的聚类信息处理和组间差异表示等[4]。由于最小生成树对全脑网络的过度简化，使用局部可量化指标表达信息较少，进而导致组间差异表达能力降低。研究表明，仅采用最小生成树网络特征来进行分类研究，较其他网络特征而言，其分类准确率明显偏低[5]。

本文提出了一种在局部差异网络的基础上构建最小生成树功能脑网络进行特征提取的新方法，可最大程度实现组间差异表征，且能提供更有效的分类特征，以服务分类研究。该方法通过对抑郁组和正常组的脑网络进行组间基于网络的统计(network-based statistic，NBS)，寻找组间功能连接强度有明显区别的连接以连接涉及的大脑区域；再分别构建以每个脑区及其差异连接涉及的脑区为节点的功能连接子网；在此基础上，对每个子网进行最小生成树功能连接网络的构建，进行进一步的分析研究。

本文研究框架如图1所示。首先采集数据并进行预处理，其次利用Pearson相关构建低序功能连接网络，在此基础上利用NBS找到差异连接，并分别构建局部差异网络，接下来采用Kruskal算法构建局部差异最小生成树脑网络，进行局部指标的计算并提取特征，最后构建分类器验证分类结果。

图1 本文的研究框架Fig.1 Research framework of this paper

1 数据采集与预处理

1.1 被试人员

本研究的受试者由山西医科大学第一医院招募，共70名。在后续的实验数据处理过程中发现，有4例受试者的数据不符合要求，因此，本文只对其余66名受试者的数据进行研究。其中，38名是首次发病且未服用任何药物的重性抑郁症患者；另外28名是健康志愿者，作为本实验的对照组。在进行扫描前，研究人员与所有被扫描者通过书面方式达成一致意见(抑郁症患者与其家属签订协议，对照组志愿者与本人签订协议)。患病严重程度由24项Hamilton Rating Scale for Depression(HAMD)以及Clinical Global Impression of Severity(CGI-S)来表征。具体如表1所示，其中，数据范围是最小值-最大值(平均值±标准差)。其中，a表示双样本双尾T检验，b为双尾皮尔逊卡方检验。

表1 被试人员基本信息Table 1 Demographics and clinical characteristics of the subjects

1.2 数据采集与预处理

本实验使用的数据均是由山西医科大学第一医院的放射科专业磁共振医师采集。使用的是德国西门子磁共振设备(siemens trio 3-tesla scanner,siemens,erlangen,germany)。受试者在扫描期间，研究人员用海绵固定其头部，以防止来回晃动。受试者需保持身心的放松，闭目，处于一种不去想特定的事情但又不能进入睡眠的状态。扫描参数设置为：回波时间(Echo Time)30 ms，射频重复时间(Repetition Time)2 000 ms，层厚(Slice Thickness)4 mm，层间间隔0，视野范围(Field of View)192 mm×192 mm，翻转角90°，存储矩阵64×64.实验共采集了248个时间点。此外，由于磁化的稳定性，磁共振扫描数据的时间序列的前十项被抛弃。

本次研究的受试者磁共振成像数据利用SPM8(statistical parametric mapping)软件做预处理。在对数据整体做头动校正和时间片校正时，抑郁组(major depressive disorder，MDD)有2例受试者的数据、正常组(normal control，NC)有2例受试者的数据因为其头部运动大于3 mm或者其头部转动幅度大于3°，在校正过程中被舍弃。本文的66例实验数据不包含此4例数据。通过12维仿射变换，使校正后的图像得到进一步优化，然后将其标准化到3 mm体素的MNI(montreal neurological institute)标准模板上。最后，用线性下降和0.01～0.10 Hz带通滤波缩小高频生物干扰以及低频漂移对实验的影响。

2 网络构建与属性计算

2.1 节点和连接的定义

本文利用MNI机构的标准脑解剖图AAL(anatomical automatic labeling)模板将大脑划分成90个脑部区域(左右半球对称，各包含45个区域)。脑网络中的每个节点用所划分的大脑区域来表示[6]，节点值是每个脑区中全部体素的算术平均值。

在本研究中，Pearson相关系数代表脑网络的边。扫描过程中会产生一定的头部移动的伪差异，脑脊液还有白质信号的噪声，通过用多元线性回归的方法减小其对实验的影响。通过Pearson相关，对残差进行计算，就得到了任两节点间平均时间序列的相关系数γij，生成一个90×90的相关矩阵。数学定义为：

(1)

2.2 局部差异网络的构建

本文使用NBS方法识别MDD和NC之间功能连接强度存在显著组间差异的连接及涉及到的脑区[7]。主要的步骤如下：

1) 对Pearson相关矩阵应用Fisher r-z转换，求Z-score矩阵，提升相关性系数的正态性；

2) 采用多元线性回归方法对Z-score进行年龄和性别混杂因子的校正；

3) 在MDD和NC的Z-score矩阵之间进行NBS.

本文将模板中的90个脑区分别作为90个ROI(Region of Interest)，以每个ROI及其与之相关的差异连接涉及到的脑区为节点，构建功能连接子网，局部差异网络构建完成。

2.3 局部差异最小生成树构建

神经精神脑疾病研究中，Kruskal算法常用于脑网络最小生成树的构建[8-9]。本文使用Kruskal算法[10]进行局部差异最小生成树脑网络LDN-MST(minimum spanning tree based on local difference network)的构建。详细过程：第一步，将皮尔逊相关矩阵里所有的边按照权重值大小从大到小依次排列；第二步，连续向网络加入权重最大的边。其中，如果加入边后有环出现，那么就把加入的这个边丢弃；当所有的边都被放入到网络中时， Kruskal算法构建过程就结束了。最后，MDD组和NC组的LDN-MST构建完成。

2.4 MST属性计算

构建完成局部差异最小生成树脑网络后，要进行网络属性值计算和网络拓扑结构的对比。经过相关文献的查阅和大量的试验分析，介数中心度、度、离心率，以及树层次、叶子分数、直径和平均离心率分别是本文最后选用的3个局部属性和4个全局属性。

此外，为了判断混杂因子(年龄、教育程度、性别)对最小生成树4个全局属性的影响，本文通过多元线性回归方法进行分析。分析表明，混杂因子与局部差异最小生成树脑网络的平均离心率、树层次、叶子分数、直径等4个属性之间没有显著相关性。

3 特征选择和分类

3.1 统计分析

为了对MDD组和NC组局部差异最小生成树网络的3个局部指标进行统计分析，本次研究采用Kolmogorov-Smirnov非参数置换检验检测组间存在显著差异的大脑区域。同时，利用该检验对网络的直径、叶子分数、树层次和平均离心率4个全局属性进行分析，并对NC组和MDD组的网络拓扑结构进行比较。

本文使用Benjamini-Hochberg假阳性率法(q=0.05)对结果进行校正。FDR(False Discovery Rate)方法适用于样本量相对小的比较结果的校正，因为总Ⅰ型出错率可以较好地在多重比较中被控制。

3.2 分类

本文利用机器学习构造分类器。被人们广为使用的分类器多种多样，如决策树、支持向量机等，根据不同的情境选择合适的分类算法。支持向量机(support vector machine，SVM)这一分类算法比较适合且经常被用于对数据量小的样本的处理[11]。根据功能磁共振数据的自有特征，本文选用SVM进行分类研究，通过径向基核函数进行分类模型的构建，使用留一交叉验证来度量分类器的泛化能力。具体而言，假设总样本数为N，每一次都从总样本中拿出其中的一个样本作为测试集，那么除此之外的N-1个样本就成为训练集。按这样的规则，确保每个样本都可以做一次测试样本，这样将会产生N个分类器并通过计算得出N个测试结果，最后N个结果取平均值来衡量分类器模型的性能。

本文用正确率、灵敏度、特异度和AUC值这4个最常被使用的度量指标作为标准。其中，正确率是样本中被正确分类的样本数除以样本总数所得的商；敏感度体现的是正例被正确分为正例的样本数与全部正例数之比，可以衡量出分类器对正例的辨别能力；特异度体现的是负例被正确分类为负例的样本数与全部负例数之比，能够衡量出分类器对负例的辨别能力。

4 结果

4.1 LDN-MST网络

如图2所示，通过基于网络的统计我们发现了MDD与NC存在有明显组间差异的连接(P<0.05，FDR校正)。具体的，有28个ROI与其余脑区间的连接无显著组差异；有26个ROI与其余脑区间存在显著组间差异的连接数为1条；有11个ROI存在显著组间差异的连接数为2条；其余25个ROI存在显著组间差异的连接数大于等于3条。表2列举了差异连接条数大于等于3的脑区及其具体差异连接条数。

图2 差异连接及脑区Fig.2 Differential connection and brain area

本文将差异连接条数大于等于3的ROI和其差异连接涉及脑区共同作为局部差异网络的节点来构建局部差异网络，并且在其基础上构建了最小生成树功能连接子网络，最后一共构建了25个功能连接子网。

4.2 全局属性

本文通过Kolmogorov-Smirnov非参数置换检验，对平均离心率、树层次、叶子分数、直径4个指标进行了分析，辨别NC组与MDD局部差异最小生成脑网络的拓扑差异(P<0.05，FDR校正)。其中，25个子网中只有6个ROI形成的子网其指标有显著性差异，且差异只出现在平均离心率和直径这两个指标上。抑郁组的平均离心率和直径均比对照组低。

表2 脑区及其差异连接数量Table 2 Brain regions and the number of different connections

全局属性检验结果如图3所示。图中为，以(a)左侧中央前回、(b)左侧中央沟盖、(c)右侧楔叶、(d)左侧豆状壳核、(e)左侧颞上回以及(f)左侧内侧和旁扣带脑回为ROI网络的直径与平均离心率的结果。其中，每幅图中左侧为直径的结果，右侧为平均离心率的结果；图(f)中，左侧内侧和旁扣带脑回为ROI的局部差异最小生成树网络中，只有在平均离心率这一指标上存在明显组差异。

4.3 局部属性

对每个受试者，计算了局部差异最小生成树脑网络的离心率、介数、度，并且使用Kolmogorov-Smirnov检验来分析计算的指标。选择MDD组与NC组之间有特别明显差异的脑区作为分类特征进行之后的实验。

如图4所示，这些大脑区域的局部属性都有明显的组差异(P<0.05，FDR校正)。

4.4 分类

本文使用留一交叉验证法(leave-one-out cross validation，LOOCV)度量分类器模型的泛化能力。评价分类器性能的指标很多，本文用准确性(Accu-racy)、灵敏度(Sensitivity)、特异性(Specificity)以及AUC(area under curve)这几个最常被使用的指标作为度量标准。

图3 全局属性结果Fig.3 Results of global properties

图4 局部属性结果Fig.4 Results of local properties

其中，选取作为分类的脑区特征是MDD和NC每组子网中的离心率、度、介数中心度3个局部属性中组间差异非常明显的总共24个特征，用来进行之后的分类，具体见表3所示。

在表4中，SACCHET et al[12]和ERGUZEL et al[13]通过全局指标进行分类；GUO et al[14]采用的分类特征是3个局部指标：节点效率、度、介数；WEE et al[15]的分类特征是局部聚类系数；GUO et al[5]在全脑构建最小生成树网络时，利用局部指标进行分类。前面所提到的分类是按脑区的相关信息进行的，多个脑区之间的拓扑信息受到损失。本文所提方法能保证一定的组间差异表征能力，能获得更多的组间信息，弥补在全脑上用最小生成树构建网络的不足，特征的有效性和分类的准确率得到进一步的提升，最高分类准确率可达83.3%.