静息态功能脑连接的空间动态分析及分类研究

高 晋，赵云芃，Godfred Kim Mensah，李欣芸，刘志芬，陈俊杰，郭 浩

1.太原理工大学 信息与计算机学院，太原030024

2.太原理工大学 艺术学院，山西 晋中030600

3.山西医科大学第一医院 精神卫生科，太原030000

静息态功能脑网络分析[1]使用基于血液氧合水平（Blood Oxygenation Level Dependent，BOLD）的功能磁共振成像技术，是研究大脑功能的一种重要方法[2]。研究表明，神经精神疾病患者的临床表现与其大脑功能网络连接异常是相关的[3]。脑网络与机器学习相结合的方法已经被广泛应用到脑疾病的诊断[4]中，如精神分裂症[5]、阿尔兹海默症[6]、癫痫症[7]等。因此，静息态功能脑网络分析方法在脑疾病的分析和诊断中非常重要。

在传统的脑网络分析中，隐含的假设是大脑功能连接在整个静息态功能磁共振扫描过程中是恒定不变的[8]。然而，无论是在经验上还是通过实验都证明了大脑功能连接随时间推移而发生动态变化[9]。Wee等使用滑动窗口的方法构建了静息态时间动态网络，并应用于早期轻度认知障碍病人的识别中[10]。

此外，近期的研究发现除了时间动态性外，脑网络存在着显著的空间动态性。目前，领域内关于空间动态的定义主要有两个方面：一是以脑区内体素BOLD信号值随体素空间位置变化的动态变化来表示脑区的空间动态性[11-12]；二是以脑区与脑区边缘体素之间的连接随边缘体素位置变化的动态变化[13-14]。同时，空间动态在提取更丰富的特征及提高分类准确率上也取得了显著的效果。Gopal等的实验结果证明精神分裂症病人与正常对照组相比，颞上回与颞中回脑区的空间动态性增大[15]。接标等提取脑区的空间动态特征，用于AD/MCI分类，发现空间动态可以提高AD/MCI的分类性能[16]。

但是，上述研究的研究对象都是脑网络的节点（成分或脑区）的空间动态，并没有研究针对脑区空间动态变化所引起的功能连接变化开展工作。同时，无论是时间损耗上还是空间存储上，每个脑区成千上万的体素数在计算脑连接空间动态时都会使得算法复杂性很高。因此，如何表示脑连接的空间动态变化来更加有效地提取特征是一个具有挑战性的问题。基于此，通过高维模板将每个脑区内体素整合成多个子脑，以两脑区子脑之间功能连接随空间位置变化的动态变化来表示两个脑区之间功能连接的空间动态。

本文提出了一种基于脑连接空间动态的分类方法。具体地说，本文以抑郁症为疾病基础，采用AAL模板进行脑区分割，计算功能连接矩阵，计算NBS（Network-Based Statistic）差异连接，提取差异连接的空间动态特征，选取KS 检验p ＜0.05 的特征作为分类特征进行分类器模型构建。结果表明：与传统的分类特征相比，研究中提取的差异连接空间动态特征分类性能更高。

1 数据

1.1 被试

本研究在实验之前招募了包括38名首发未用药的抑郁症（MDD）患者和28 名健康被试共66 名被试。患者都为山西医科大学第一医院精神卫生科所确诊的中国籍汉族抑郁症患者。本研究在山西省医学委员会的建议之下，按照世界医学协会赫尔辛基宣言与所有被试签署了书面协议。所有被试基本信息统计如表1所示。

表1 被试基本信息统计

1.2 数据采集与预处理

数据的采集工作在山西医科大学第一医院由熟悉磁共振操作的放射科医生来完成。所有被试都在3T超导磁共振扫描仪（Siemens Trio 3-Tesla scanner，Siemens，Erlangen，Germany）下进行静息状态下的功能磁共振成像（fMRI）。在扫描过程中，被试头部被海绵固定以防止产生头动位移，同时要求参与者放松，闭上眼睛，不去想特定的事情但不能入睡。每个扫描包括248 个连续的EPI 功能图像，扫描参数设置如下：射频重复时间为2 s，存储矩阵为64×64，回波时间为30 ms，层厚与层间间隔分别为4 mm 和0 mm，成像视野为192×192 mm。由于初始磁共振信号的不稳定性以及被试对环境的自适应性，前十个功能图像的时间序列被丢弃。

图像预处理使用的是SPM8工具包（http：//www.fil.ion.ucl.ac.uk/spm）。首先对数据集进行时间片校正和头动校正。剔除头动大于3 mm 或转动大于3°的被试数据，最终获得了66例被试数据。然后，图像进行12维度的优化仿射变换，将其标准化到3 mm体素的MNI标准空间中。最后，进行线性降维和低频滤波处理。

2 方法框架

针对静息态脑连接的空间动态分类方法，图1给出了该方法的框架，包括以下步骤：第一步，数据采集与预处理；第二步，脑网络构建，计算每个被试的功能脑连接矩阵，使用NBS 方法获取组间差异连接；第三步，特征提取及选择，这里提取了脑连接空间动态特征，使用KS检验进行特征选择；最后，SVM 训练与预测，使用10 折交叉验证进行分类器模型构建。

图1 本方法的框架

2.1 脑区分割

本研究在数据预处理后，根据AAL（Anatomical Automatic Labeling）模板将大脑划分为90个（左右半脑各45 个）脑区。通过计算每个脑区内所有体素在不同时间点的BOLD信号值的算术平均值来获得90个脑区的平均时间序列。

为了进一步表示脑连接的空间动态，本研究根据基于AAL模板的独立脑区进一步划分得到的Parc1501模板将大脑划分为1 501个脑区[17]。然后，提取1 501个脑区平均时间序列。这里的Parc1501 模板的每一个脑区只且一定包含于AAL模板的一个脑区中。例如，AAL90模板下的左侧中央前回在Parc1501模板下被划分为32个规模更小的子脑。两模板的节点映射如图2。

图2 两模板节点映射图

2.2 脑网络构建

在获取到AAL 模板90 个脑区的时间序列后，本研究以90个脑区作为节点，以节点之间的功能连接（Pearson相关系数）构建功能连接网络。对于每个被试，本研究获得一个90×90的加权全连接网络。第i 个脑区与第个j 脑区之间的Pearson相关系数的数学定义如下：

其中，xt、yt(t=1,2,…,T)表示第i、j 个脑区的时间序列在时间点t 处的值，xˉ、yˉ表示第i、j 个脑区的时间序列的均值，T 为时间序列的长度。

为获取抑郁组和健康对照组之间的差异连接，本文使用了NBS的方法。该方法为每条连接计算了统计检验值PNBS，本文选取PNBS小于0.05的连接作为抑郁组和健康对照组之间的差异连接。

2.3 脑连接的空间动态

如图3 所示，AAL 模板下的第i 个脑区A 与第j 个脑区B 在Parc1501 模板下的M 个子脑与N 个子脑的时间序列分别为x={x1,x2,…,xm,…,xM} 以及y={y1,y2,…,yn,…,yN}。对于第i 个脑区A与第j 个脑区B之间的脑连接ρi,j，分别计算时间序列xm(m=1,2,…,M)与yn(n=1,2,…,N)之间的皮尔逊相关系数，得到M×N条脑连接的连接矩阵ρi,j(m,n)∈[-1,1](m=1,2,…,M;n=1,2,…,N)。本研究以这M×N 个连接值的动态变化来表示脑连接ρi,j的空间动态。研究中将这M×N个值降维为组距为0.1，组数为20的直方图，以直方图的高度变化来表示脑连接的空间动态变化。研究中将直方图的高度作为分类器的分类特征，K 条差异连接的可用来分类特征总数为20K 。

图3 脑连接的空间动态图示

本研究利用Kolmogorov-Smirnov 非参数置换检验对抑郁组和正常组的脑连接的空间特征进行显著性检验，以得到组间差异显著的特征。此后，利用Benjamini-Hochberg 假阳性率法(q=0.05)对结果进行校正。该方法能够在多重比较中较好地控制总I 型错误率，更适合小样本比较结果的校正[18]。

2.4 特征选择与SVM分类

为了找到特征的最优子集，避免过度拟合、提升模型性能、更快地训练分类器，需要在分类前进行特征选择。本次研究选择统计显著性P 值作为分类特征选择方法（P ＜0.05,FDR校验）。

本文采用支持向量机（SVM），这是相似的研究中常用的方法[19-20]。特别是其对小样本数据，具有良好的分类效果[21]。本研究中，选择SVM 作为分类器是基于MATLAB 的LIBSVM 工具包（http：//www.csie.ntu.edu.tw/～cjlin/libsvm）进行分类。SVM核函数为径向基核函数（Radial Basis Function，RBF），其中，惩罚因子C=1，核参数g=0.1。

为了研究脑连接的空间动态对脑疾病分类器性能的影响。本文分别以脑网络的K 条差异连接值、差异连接的空间动态特征作为分类特征进行分类器的构建。本文使用10折交叉验证（10-fold cross validation）[22]的方法来评估分类器的泛化性能。具体的过程是将所有的被试随机分成10等份，逐一将其中的1等份作为测试集，剩余的9 等份是训练集，最后对10 次结果的均值作为对分类器性能评估。同时，为了得到更精确的结果，本实验进行100 次10 折交叉验证，最后对100 次的结果求均值得到最终的结果。

2.5 分类特征评估

为了评估所选特征与分类器的关联性，采用了ReliefF算法[23]。ReliefF是一种有效的特征选择算法，根据各个特征和类别的相关性赋予特征不同权重。特征权重越大，表示该特征分类能力越强，反正，表示分类能力越弱。本文采用ReliefF算法来量化两组分类特征的重要程度，并用K-S检验进行显著性分析。

3 实验结果及分析

3.1 差异连接结果及分析

每个被试的脑连接矩阵为90×90，因此每个被试的脑网络中有4 005 条脑连接。构建完脑网络后，本文使用NBS 方法对66 个被试的脑网络进行了统计分析。NBS 会为每一条脑连接计算一个PNBS值，PNBS＜0.05的脑连接为具有组间差异的脑连接。本文共获得87条PNBS＜0.05 的脑连接（图4）。

图4 87条脑差异连接

与对照组相比，MDD 患者的一些功能连接增强（如，左侧丘脑与左回直肌之间的功能连接，左侧丘脑与右侧眶内额上回之间的功能连接等），一些功能连接减弱（左侧杏仁核和左侧内侧和旁扣带脑回，左侧颞中回和左侧中央沟盖等）。本文按照PNBS筛选出前10条差异连接（表2）。

表2 NBS差异连接top10

这些脑差异连接在以往的研究中也有类似的结果。Zhou等[24]采用静息状态FC分析方法，对18例MDD患者和20例健康对照者来识别内在组织的脑区和研究与单个脑区相关的脑连接。结果表明，抑郁症患者的前扣带回与双侧前额叶之间的脑连接减弱。Greicius等[25]利用独立成分分析方法，发现抑郁症患者额叶眶部与丘脑之间的功能连接强度显著增加。Ichikawa等[26]在静息态fMRI数据的基础上，开发了分类器来区分MDD患者和健康对照组。研究中指出，在补充运动区与岛盖部额下回之间的功能连接为分类中所占权重最大的十二条脑连接之一。

另外，与这些脑差异连接相应的脑区也被认为是抑郁症诊断中比较重要的脑区。Qin等[27]应用复杂网络方法，探讨了脑解剖网络的拓扑结构，并指出右前额包括右额中回与右回直肌的中心度与抑郁症的病程呈负相关。郭浩等[28]利用基于Elastic-net 方法构建了超网络，并利用统计计算方法得出：右侧补充运动区为十二个异常脑区之一。另外，其他脑区包括（部分边缘系统（左侧杏仁核，内侧和旁扣带脑回，左侧前扣带和旁扣带脑回），部分颞叶（颞上回），部分额叶（岛盖部额下回，眶内额上回），部分纹状体（左侧豆状壳核，左侧豆状苍白球））均为边缘系统-皮层-纹状体-苍白球-丘脑神经环路（LCSPT），而LCSPT 环路被广泛认为是抑郁症的主要病理环路[29]。

3.2 空间动态特征结果

在获取到抑郁组与正常对照组之间的差异连接后，对脑连接进行了空间动态分析，并分别提取了抑郁组和对照组脑连接的空间动态特征。然后利用KS检验的P值来选取最具有判别性的空间动态特征。本研究将P值阈值设置为0.05，得到最具判别性特征124个。

另外，为了更好地研究抑郁组与正常组之间的异常连接，本研究按照脑连接包含最具判别性空间动态特征的数量进行排序。选取空间动态特征数量较多的脑连接进行分析。与NBS 差异连接相同的是，双侧岛盖部额下回之间的脑连接的空间动态也为最具判别性特征的脑连接。另外，左侧豆状苍白球和左侧缘上回之间的脑连接、右侧内侧和旁扣带脑回和左侧岛盖部额下回之间的脑连接、双侧顶上回和左侧补充运动区之间的脑连接、左侧杏仁核和左侧内侧和旁扣带脑回之间脑连接的空间动态特征在抑郁组和正常组之间的差异更显著。

3.3 分类结果分析

传统的基于静息态功能脑网络的脑疾病分类方法在特征提取时，提取一些脑网络的全局属性、局部属性等作为分类特征。为了证明本研究中基于脑连接空间动态的分类方法具有较好的分类性能，说明研究中提取的脑连接的空间动态特征的优越性，本文对比了不同特征提取方法的分类准确率、敏感性、特异性和AUC。正如表3所示，与已有的提取脑网络静态特征的脑疾病分类方法相比，脑连接的空间动态分类方法分类性能都较高。另外与已有的基于脑区的空间动态分类方法相比，本研究中的方法也取得更高的分类准确率。

表3 不同方法的分类性能

除此之外，上述不同方法的研究使用的实验数据不同，本研究方法得到较好结果可能是数据集的影响。为了更准确地评估脑连接空间动态对分类的影响，本研究同时构建了基于差异连接的分类器。分类结果表明：基于脑连接的空间动态特征的分类方法的分类准确率提高了5.2%。

最后，如图5，脑连接的空间动态特征的reliefF的平均值较大，分类能力更强。其潜在原因在于，与传统的静态功能连接特征相比，脑连接的空间动态分类方法在特征提取步骤有效地捕获到脑连接的空间动态信息。

图5 两组特征的reliefF的平均值

综上所述表明，本研究的静息态脑连接的空间动态分类方法可以有效地用于正常人与抑郁症病人的分类中。

4 结束语

静息态功能脑网络的研究已经为神经性疾病的检测与诊断带来了极大的帮助，但现有的研究主要集中在大脑之间功能连接的静态特征，忽略了脑连接的空间动态信息。为解决这一问题，本研究提出了一种静息态功能脑连接的分类方法。具体地说，首先进行数据采集与预处理，然后构建脑网络及获取NBS差异连接，提取脑连接的空间动态特征，特征选择，最后SVM分类模型构建。分类结果表明，静息态功能脑连接的空间动态分类方法取得了可观的分类效果。但是，该方法存在一定的局限性，其一，本研究只采用AAL 一种分割模板，没有研究此方法在其他模板中是否会得到同样的结果，未来工作中将会使用不同模板评估本文所提出的方法；另外，该方法中只针对差异脑连接的空间动态进行了分析，虽然降低了算法运行时间，但可能忽略了其他潜在的脑连接的空间动态信息，因此怎样在不增加算法消耗的情况下，进一步改进该方法是今后研究的方向。