动态脑功能网络Rich-club时空观测方法

盛景业，王 彬，2，薛 洁，淡杨超，刘 畅

（1.昆明理工大学信息工程与自动化学院，昆明，650500；2.昆明理工大学云南省人工智能重点实验室，昆明，650500；3.云南省公安厅禁毒局，昆明，650228）

引 言

大脑是人体最复杂的动态信息处理系统，利用各种技术方法对大脑功能进行探究是目前的研究热点之一。其中功能磁共振成像（Functional magnetic resonance imaging,fMRI）成像是一种非侵入式神经影像学技术[1]。其原理是通过磁共振成像技术采集fMRI-BOLD信号，并基于BOLD信号观测大脑血氧水平的变化，以此来探究大脑的活动情况[2]。在以往的相关脑网络研究中[3-4]，有些研究者认为时间分辨率为无穷大，即把fMRI-BOLD信号的采集时间段内的大脑连接取平均后合成为一个状态，但近年来大量研究表明，大脑活动是一个随时间变化的动态过程，脑功能网络及BOLD信号变化可以反映出脑功能网络的动态属性[5]。基于fMRI信号技术重构的脑网络是以大脑的不同脑区作为网络节点，两两脑区之间的某个属性作为网络的边构建而成，它是一个复杂网络系统，其中存在着一类核心脑区节点，这些脑区的数量较少，但是与其他脑区的联系较为紧密，并且影响力较大，往往在整个脑网络中占据了举足轻重的地位。这些脑区的集合一般被称为脑区的Rich-club结构[6]。脑网络中重要脑区的识别以及Rich-club属性的评价是脑科学研究中的重要方向，研究者期望根据复杂脑网络中的核心节点和Richclub结构来进行大脑的运行机理研究和脑疾病的辅助诊断。Guusje等[7]采用基于节点度的评价方法对精神分裂症及其亲属的核心节点和Rich-club属性展开研究，发现精神分裂症患者相较于正常人来说其脑网络中节点度较大的脑区受损更为严重，同时病人亲属的脑区属性与正常人相比较也有所下降。Geetharamani等[8]采用节点度特征向量中心性、介数中心性等方法对大脑的不同脑区进行重要程度评价，发现丘脑、岛叶和海马等脑区重要程度较高。Jin等[9]采用不同的先验模板进行脑区分割，并选取不同的重要程度评价方法来标定重要脑区，发现样本在使用相同的模板分割脑区后，无论采用哪种重要程度评价方法所得重要脑区大体相同；与此相反，若采用不同模板划分脑区，即使采用相同评价方法，所得重要脑区也明显不同。Martijn等[10]将静息态脑网络进行模块化处理，再将各模块内重要脑区进行标注。然而此方法只能判断对模块内部来说较为重要的脑区，忽略了模块之间连接脑区对于整体脑网络联通性的影响。Weng等[11]基于重要度贡献矩阵判定人脑网络核心节点，分别判定了无权脑网络核心节点及加权脑网络核心节点，并证明该方法的合理性。

上述研究主要致力于通过判定静态脑结构网络的核心节点，将核心节点的属性作为重要指标来进行人脑功能研究或脑疾病判断。然而越来越多的研究围绕着基于fMRI重构的动态脑功能网络而展开，这是因为动态脑网络不但包含脑区的空间属性，还包含了时间信息。目前大多数动态脑网络的研究着重于全脑网络的观测方法。Wee等[12]通过对动态脑网络的稀疏特性分类来对轻度认知障碍做出早期诊断；O″neill等[13]将静息态网络利用固定小时间窗口重组，得到动态脑功能网络并对比静态和动态网络的区别；杜宇慧等[14]提出一种基于独立成分分析法进行动态脑功能网络分析的方法，比较正常被试和精神分裂患者在功能连接状态的差异。Lindquist等[15]采用动态条件相关模型(Dynamic conditional correlation，DCC)发现脑网络中不同脑区功能连接并非静止不变，由此确定脑网络是一个动态变化的过程。Li等[16]采用K-means算法对fMRI数据降维并提取特征；由于K-means算法不能完全映射出样本点在未降维前真实分布情况，且异常点对均值影响较大，因此刘畅等[17]提出一种基于瞬时转变率模型的脑网络状态降维及观测算法对脑网络状态转换趋势进行观测。

由于基于fMRI重构的动态脑功能网络是随时间变化的动态连接网络，其中每一个时刻脑网络的核心节点是不同的，因此，Rich-club的组成也会随时间而变化。由于Rich-club是所有脑区中最重要的节点集合，通过观测Rich-club组成随时间的变化，得到脑网络空间结构与时间变化的联系，就能将脑区的空间属性与时间属性相观测，这对于研究脑网络的状态表达和动态属性具有重要意义,但是目前此类研究较少。本文以基于fMRI-BOLD的动态脑网络作为研究对象，以人脑Rich-club的辨识技术为基础，以观测重要脑区组成随时间的变化为目标，提出了一种动态脑网络Rich-club时空观测模型的构建方法。

1 Rich-club时空观测模型构建方法

动态脑网络的Rich-club时空观测模型构建方法原理如图1所示。首先，将原始的fMRI-BOLD数据进行预处理后，按照AAL（Anatomical automatic labeling）模板[18]提取90脑区在若干时间点上的BOLD信号；采用滑动窗口技术和Person[19]计算法计算不同脑区之间随时间变化的相关系数矩阵，得到信号采集时间上的多个脑功能网络。然后根据核心节点判定方法对每一个脑功能网络中的核心节点进行标记再组成Rich-club集合，通过对样本内所有时间点的Rich-club进行聚类而将空间与时间观测在一起，并构建观测时间区间上的Rich-club重要性占比指标，用以定量评价动态脑功能网络的重要性。

图1 Rich-club时空观测模型构建方法Fig.1 Construction method of Rich-club spatiotemporal observation model

1.1 瞬时脑功能网络核心脑区定位及Rich-club构建

瞬时脑功能网络核心脑区定位是指对每个时间点所重构的脑功能网络中的重要节点进行辨识。由于动态脑功能网络是由每一个采样时间点上的相关系数组成的加权网络，因此应选用合适的加权核心节点判定方法[20-22]。由于本文的研究重点是动态脑功能网络中Rich-club组成随时间变化的情况，只需确定Rich-club中核心节点的组成情况，暂时不需要分析核心节点之间的连接关系以及节点在网络中的位置信息等因素。因此本文采用加权节点度法作为核心节点的判定方法，所使用的网络参数为加权网络中节点度，即节点与其他节点之间连接权重之和。该方法还具有最能体现网络中节点连接重要性、计算量小以及复杂度最低等优点。

加权节点度法的计算公式为

式中：di为脑区i的节点度，wi,j为脑区i与脑区j连边的加权权重值。由式（1）可知节点度di即为该脑区与脑功能网络中其他各脑区的连接权重总和。di值越大，证明脑区i在该网络中重要性越高。由于本文采用AAL模板[18]将脑网络划分为90个脑区节点，因此计算脑区i节点度时，需要将脑区i与其他89个节点之间的连接权重累计相加，在计算时视脑区与自己本身之间的连接权重为1，同样计入节点度中。脑功能网络判定核心节点过程如图2所示。

按照式（1）计算出网络中各脑区在某一特定时间上的节点度大小后，选取节点度最大的前10个节点组成该时间点的上Rich-club[23]。

图2 脑功能网络判定核心节点过程Fig.2 Brain functional network to determine the process of core nodes

根据脑功能网络的时间先后顺序，可以得到所有时间上对应的Rich-club集合，由于大脑的状态会随着时间的推移而产生变化，因此脑区之间的连接强度也会发生变化。同时，这将引起核心脑区组成成分发生变化，为了观测脑功能网络的核心脑区随时间变化关系，以Rich-club集合为对象展开进一步的研究，即构建脑功能网络的时空观测矩阵，以观测不同时间点上脑功能网络Rich-club的相似性，从而探究脑功能网络随时间的演化模式。

1.2 Rich-club时空观测模型

与静态脑网络的研究不同，由于脑功能网络的Rich-club结构随时间发生变化，其组成的核心节点也将发生变化，但总体说来，其变化模式应该具有一定的规律性。因此以1.1节得到的Rich-club为对象，对其随时间的演化模式展开研究，并构建时空观测模型。

时空观测模型的构建过程如图3所示。首先将1.1节中所得每个时间点上的Rich-club集合转化为一个10维Rich-club向量，并计算各时间点的Rich-club脑区相似性，然后采用聚类算法对一个被试样本在所有时间点的脑功能网络的Rich-club进行状态分析，观察其状态变化与时间之间的变化模式。

图3 可视化分析步骤Fig.3 Visual analysis steps

本文采用Dice系数来衡量各时间点的Rich-club脑区相似性。Dice系数可用于衡量高维空间内两向量的相似性，其原理为计算两时间点的Rich-club向量中相同的核心脑区数量，再与每个时间点上的脑区Rich-club向量长度进行对比，所得的比值即代表了两个时间点上脑区Rich-club的相似程度，Dice系数取值范围在0～1之内，Dice系数的计算公式为

式中：comm(ti,tj)为时间点ti向量和时间点tj向量的相同核心脑区的个数；leng(ti)和leng(tj)为ti，tj向量的长度。

计算各时间点Rich-club脑区的Dice系数后，即可构建出一个N×N的脑网络Rich-club时空观测矩阵，其中矩阵中各数值所表示的含义为某两个时间点之间Rich-club的相似程度。该矩阵能够表达在不同时间点上的重要脑区组成的相关程度和变化情况，为了探究人脑在不同时间点上Rich-club的变化情况，将Dice系数作为不同时间点之间的空间相似度指标，利用DBSCAN(Density-based spatial clustering of applications with noise)聚类算法将所有时间点的Rich-club向量进行聚类。

DBSCAN算法是一种基于密度的空间数据聚类方法[24]，该算法可以将样本空间内具有足够密度的部分区域划分为一簇，并且可以将样本数据中不满足聚类条件的噪音时间点筛选出来，这样能够避免脑网络采集过程中的噪音点对聚类效果的影响。通过DBSCAN算法可以将每一个被试样本在多个时间点上的Rich-club集合聚为n类，其中纵坐标为-1的数据表示的是噪音点，其余纵坐标数值表示的是不同状态类型。该聚类结果不但显示出脑功能网络随时间变化过程中其状态也发生了改变，也显示出Rich-club在相邻时间内的一致性，能够观测到状态的转换，这些结果可以直观地表达出脑功能网络随时间变化状态逐渐发生变化的过程，完整地表达了在数据采集期间内脑网络中关键脑区的演化模式。

1.3 动态Rich-club脑区重要性评价模型及指标

目前的脑区重要性评价指标一般只考虑某一种因素[25]，上文中所述时空观测矩阵和Rich-club的演变模式都是基于多个时间点上连续变化的分析，由于观测了时间与空间的信息，每一个实验样本都将基于多个时间点而产生多个脑功能网络，在不同时间层上组成Rich-club的脑区也将不断地发生变化，因此单一的脑区重要程度指标有一定的缺陷。为了定量地描述每一个被试的动态Richclub属性，本节在1.2节基础上构建一个描述全脑在整个时间区间上重要脑区的动态Rich-club评价模型。

构建模型时若只考虑各脑区节点的总和，会出现某一脑区在少数瞬时脑网络中节点度较大，而在大多时间点节点度较小的情况，会影响分类结果核心节点组成的平均值；若只考虑各脑区作为Richclub成员出现的次数，会出现几个脑区出现次数相同的情况，即无法比较这些脑区在分类中的重要程度。因此，本文将脑区权重与其作为Rich-club成员出现的频率结合，作为评价动态脑功能网络中的某个脑区重要性模型，有

式中：F为样本中总时间点个数；di为脑区在该时间点上作为核心节点的节点度大小，若脑区在该时间点上不属于核心节点，则di为0，否则di为其节点度大小；fi为脑区作为Rich-club成员出现的次数。

为了消除不同被试样本之间的个体差异性带来的影响，本文统计各样本内所有脑功能网络的总节点度之和作为相关系数，使不同样本之间的Rich-club重要程度可以相互比较，有

式中：D为样本内N个时间点上脑功能网络总节点度之和为所有重要程度大于d的脑区中总节点度之和。所得P即为样本Rich-club脑区在脑功能网络中重要性占比。

综上所述，基于时空观测矩阵的Rich-club状态聚类以及动态Rich-club脑区重要性模型分别从脑功能网络中动态Rich-club的时间相似性与动态Rich-club的重要性两个方面给出了定量的描述，为研究脑网络的状态表达和动态属性提供了基础，并为脑疾病患者与健康人在Rich-club动态特征上的比较和分析提供了可能性。

2 实验及结果分析

2.1 数据来源及预处理

本文实验采用了脑网络公开数据库International Neuroimaging Data-Sharing Initiative(INDI)[26]中的斯坦福大学所提供的40组脑网络样本。40组样本分为两类：一类为20个自闭症（Autistic spectrum disorder,ASD）患者样本，样本年龄在7.5～12.9岁之间；一类为20个健康人（Test control group,TC）对照组，样本年龄在7.8～12.4岁之间。在数据采集过程中，所有参与者保持非睡眠闭眼状态。扫描重复时间 TR（Repetition time）为 2 s，回波时间 TE（Time of echo）为 30 ms，采集矩阵为 6个 4×64，视野FOV（Field of view）为20 cm×20 cm，翻转角为80°，每个样本各采集180个时间点数据。得到rs-FMRI原始数据后，采用Python/FSL Resting State Pipeline Platform平台工具[27]进行预处理，并采用AAL模板将全脑共划分成90个脑区。处理流程包括：（1）剔除各样本前4个不稳定数据时间点；（2）时间层校正；（3）头动矫正；（4）颅骨去除；（5）空间标准化；（6）带通滤波；（7）脑区平均时间序列的粗体信号提取。

由以上所述步骤得到经预处理过的样本数据后，每个样本包含了176个时间点的数据。再利用滑动窗口技术对fMRI-BOLD信号进行分区采集，取滑动窗口长度W=30，步长L=1，得到147个时间点的观测窗口，接下来计算每个时间点上各脑区之间的BOLD信号的皮尔逊相关系数作为两个脑区的强度相关系数，即每个时间点上可以得到1个90×90的脑网络状态观测矩阵。矩阵中每个数据代表两两脑区之间的相关系数，取值范围在-1～1之间，数据绝对值越大表示两脑区相关性越强。

由于本文实验中不考虑脑区之间连通的方向性，因此将矩阵中所有数据取绝对值处理。根据现有研究[28]可知，由皮尔逊相关系数组成的脑网络相关系数矩阵为了保持较好的小世界性[29-31]以及高效性[32]，其阈值应满足在0.45～0.55之间，本文选取0.48作为阈值，即将相关系数低于0.48的两脑区视为不连通，相关系数高于0.48的两脑区连接及权重保持不变。

2.2 实验设计

以2.1节中数据处理得到的脑功能网络为对象，分别对两类40组被试样本数据展开实验分析，实验方案及步骤如图4所示。

图4 实验步骤图Fig.4 Experimental procedure diagram

依照1.2节所述方法对上述实验数据进行处理。首先确定样本中各时间点的脑功能网络的核心节点及Rich-club组成，再按照时间先后排序，将所有时间点的核心节点组成转为10维无序二值化向量，接着通过计算时间点之间的Dice系数，构建出样本的Dice系数时空观测矩阵并进行可视化分析。再利用DBSCAN聚类算法对Dice系数矩阵进行聚类，最后根据聚类结果统计Rich-club在样本中的聚类个数以及不同样本之间的对照情况。对比实验包括以下3个内容：动态脑网络Rich-club时空观测矩阵的单样本对比分析；基于Rich-club观测矩阵的单样本聚类对比分析以及动态Rich-club脑区重要性评价指标的组平均分析。

2.3 对比试验及结果分析

2.3.1 动态脑网络Rich-club时空观测矩阵的单样本对比分析

按照1.1节和1.2节中所述的方法，首先利用加权节点度法计算每个样本在各个时间点的脑区重要性，然后选取前10个重要性最高的脑区作为Rich-club成员，分别得到了147个时间点上的Rich-club组成情况，接下来计算各时间点之间的Dice系数，构建出基于Dice系数的观测矩阵。

由于脑网络构建方法的原理，不可避免地存在着样本的个体性差异，因此基于Dice系数的观测矩阵结果并不相同。但通过对两类样本的分析，可以观测到类内的相似性和类间的差异性，由于文章篇幅限制，本文采用简单随机抽样方法，每类选取3个样本，它们的观测矩阵如图5所示。图中深蓝色区域表示两时间点Dice系数较小，即相似程度不高，而黄色部分表示两时间点之间Dice系数较大，相似程度较高。由图5可以看出动态脑网络的Rich-club随时间展现出一个渐变的过程，并且在某些时间点上有明显的分段情况，可以看出以下规律：

（1）健康对照组内样本时空观测矩阵中Rich-club的状态数量普遍较少，而ASD患者组内样本时空观测矩阵中Rich-club的状态数量普遍较多；相比之下，ASD患者组状态变化更频繁。

（2）健康对照组内样本在每一个相同段内时间点分布的个数较多，而ASD患者组内样本在每一个相同段内时间点分布的个数较少，相比之下，ASD患者组状态更不稳定。

（3）健康对照组内样本在Rich-club时空矩阵内不同模块内之间间隔较为明显，而ASD患者组样本时空观测矩阵中Rich-club之间间隔不明显，相比之下ASD患者组状态的界限更加模糊不清。

图5 基于Dice系数的脑网络时空观测矩阵Fig.5 Temporal and spatial observation matrix of brain network based on Dice coefficient

2.3.2 基于Rich-club观测矩阵的单样本聚类对比分析

按照1.2节中的方法，利用基于密度聚类算法DBSCAN将动态脑网络中所有时间点的Rich-club进行聚类。而在使用DBSCAN的方法对脑网络状态进行聚类时，文中用0，1，2等数字表达聚类结果得到的不同类别，由于DBSCAN算法可以将异常噪音点筛选出来，为了将这些异常点与其它聚类结果进行区分，在进行结果可视化过程中，将异常点标注为-1，从而可以得出不同样本在整个数据采集区间上的聚类情况如表1,2和图6,7所示。

从图6和图7中可以看出，TC样本所得状态聚类结果数目较少，一般只分为两类，但是在每一个状态类内的时间点个数较多；ASD样本所得状态聚类结果数目较多，并且不同的ASD样本所得的状态类别结果不同。由此可见：（1）健康对照组内样本重要脑区的组成更加稳定，状态变化不大，在相邻时间内能够较好地保持其Rich-club结构；而ASD患者组样本重要脑区的组成不够稳定，状态变化更多，在观测时间区间上Rich-club结构的数量较多。（2）健康对照组样本聚类结果一般只有两类；而ASD患者组内样本的聚类结果类数较多，并且没有显示出统一的组数。

表1 部分健康对照组样本聚类情况Table 1 Clustering of samplesof part of healthy control group

表2 部分ASD患者样本聚类情况Table 2 Clustering of some ASD patients

图6 健康对照组两组随机抽样聚类情况Fig.6 Random sampling clustering of two groups in the healthy control group

图7 ASD患者组两组随机抽样聚类情况Fig.7 Random sampling clustering of two groups of ASD patients

2.3.3 动态Rich-club脑区重要性评价指标的组平均分析

考虑到样本个体差异性可能带来的影响，按照1.3节中所构建的脑区重要性评价指标，健康对照组与ASD患者组在各时间点上Rich-club之间连接边数与当前时刻整体脑功能网络内连接边数的比值，分别将20组患者样本以及20组健康对照组样本进行组平均后得到两组组平均数据，所得结果箱型图如图8所示。

考虑到采用组平均方法时脑网络的阈值可能带来的影响，同时考虑到脑网络的小世界性以及连通性[28]，分别从0.4～0.7逐渐改变脑功能网络中阈值的大小，并观察P值的变化，将ASD组平均数据与健康对照组组平均数据所得的动态Rich-club脑区重要性评价指标结果进行比较，其折线图如图9所示。从折线图可以得出以下结论：

（1）阈值设定大于0.45小于0.7时，在整个数据采集区间上，健康对照组的Rich-club在整个网络中所表现出的重要程度均高于ASD患者组。

图8 动态Rich-club脑区重要性评价指标的组平均值对比Fig.8 Comparison of group mean values of importance assessment indicators of dynamic Rich-club brain regions

（2）ASD组的平均动态Rich-club脑区重要性评价指标比TC组的平均动态Rich-club脑区重要性评价指标数值降低了17.70%。

（3）在Rich-club组成脑区中，系数较高的连接关系占比更大，而ASD组中系数较小的连接关系较多。

从动态Rich-club脑区重要性评价指标的组平均实验结果分析可知，相对于健康人对照组来说，由于疾病对脑功能网络的影响，在静息态实验条件下，ASD病人的脑功能网络中较为重要的脑区之间的连接关系受到了破坏，并且随时间的变化其Rich-club集合的整体动态性能均有所下降。

图9 Rich-club选取不同阈值情况下重要性指标Fig.9 Importance indexes selected by Rich-club under different threshold values

3 结束语

本文针对现有动态脑功能网络中重要脑区的变化规律这一问题展开研究，给出了一种动态脑功能网络Rich-club时空观测方法，该方法首先辨识时间采样点上的重要脑区和Rich-club集合，然后采用Dice系数对fMRI-BOLD采样周期内所有时间点上的Rich-club结构相似性进行评价再进行聚类处理，从而将脑功能网络的重要脑区集合在空间和时间融合在一起，并提出了一种全局的动态Rich-club重要性评价模型，用于定量地描述动态脑网络核心脑区的的重要性和时变特征。该模型既能体现脑区在脑网络动态变化过程中所有时间采样点上作为Rich-club成员出现的频率，又保留每一个脑区节点度对其重要性的影响。采用本文方法对脑网络公开数据库INDI上提供的ASD患者组和TC组展开了3个对比实验，尽管动态脑功能网络重构时存在着一定的个体化差异，但对于相同组内样本而言实验结果仍具有普适性，而不同组间样本之间实验结果存在一定的差异性。其中动态脑网络Rich-club时空观测矩阵的单样本对比分析实验通过可视化结果证明了TC人群比ASD患者的Rich-club结构更稳定；基于Rich-club观测矩阵的单样本聚类对比分析通过可视化结构证明了TC人群与ASD患者相比，其Richclub状态数量一致，且变化较少；动态Rich-club脑区重要性评价指标的组平均分析则通过定量分析证明了TC人群组的Rich-club在整个脑网络中重要程度高于ASD组，从而揭示了ASD患者Rich-club脑区的重要连接可能因疾病而遭到破坏的现象。上述实验结果证明了本文所述方法的普适性及有效性。

本文所给出的基于时空观测矩阵的Rich-club演变模型为研究脑功能网络的动态属性提供了一种有效表达方法，动态Rich-club脑区重要性评价指标为Rich-club的时变特性提供了一种定量的判断方法。同时本文也通过实验证明了自闭症患者在动态脑功能网络的异常表现，从而为分析对比健康人与脑疾病患者之间的差异提供了一种辅助判断方法。下一步将在本文的研究基础上进一步探究Rich-club中重要脑区模式演变的规律，同时将关注重点放在ASD人群的个体差异性上，以本文的模型为基础展开研究和分析，拟探究如何通过时空融合的方法来得到ASD患者的脑活动异常的证据和判定方法。