CN111612746B - 一种基于图论的功能脑网络中枢节点的动态检测方法 - Google Patents

Abstract

本发明公开了一种基于图论的功能脑网络中枢节点的动态检测方法。本发明在多变量中枢节点检测方法的基础上进行了改进，使得能够检测出更加可靠的符合神经科学认知活动的中枢节点。首先，利用滑动窗口的技术，把血氧信号以时间为维度平均分成几段。在每一段时间的滑动窗口内，检测出对应时间窗口内的中枢节点，从而得到一条中枢节点随着滑动窗口移动的一条变化轨迹。最后，把这条变化轨迹作为一种约束作用在多变量检测的方法上，从而能够动态的检测出更加可靠准确的中枢节点。

Description

一种基于图论的功能脑网络中枢节点的动态检测方法

技术领域

本发明涉及神经科学脑网络研究领域，具体涉及一种基于图论的功能脑网络中枢节点的动态检测方法。

背景技术

静息状态核磁共振技术(fMRI)提供了一种非侵入性的方法来测量脑血氧的变化。在静息状态下，受试者没有执行任何明确的任务,受试者会自发的产生神经活动，神经活动的波动是跟血氧浓度信号的变化相关的。因此血氧浓度信号被用来计算脑区之间的连接情况，从而构建功能脑网络。

脑网络可以分成很多模块，有些模块负责视觉，有些负责听觉等，模块化结构使我们可以更加细致地区分脑区节点的不同角色和地位。比如，有些节点在其所在的模块内非常重要,但对整个网络而言却未必重要，这类节点被称为区域性核心节点(provincialhub)，而另一些节点虽然在其自身的模块内作用有限，但它们却连接着不同的模块，维系着整个网络的连通性，这类节点被称为中枢节点(connector hub)。中枢节点连接着脑网络中不同的功能模块，由于中枢节点的高度中心性，在大脑内的联系中起着非常重要的作用，比如对信息的整合，参与到多种多样的认知活动，这些特性表明当中枢节点遭受疾病的攻击时，更容易对我们的大脑产生影响。最近，在神经影像领域达成一种共识，即在整个扫描时间内，即使在无任务的环境中，大脑的功能网络也会发生变化。许多研究表明动态模式与某些脑神经疾病更相关。因此，如果我们能更加准确的动态检测中枢节点，可以帮助我们更好的分析和理解这些疾病的病理机制，也可以用来辅助疾病的早期诊断和治疗。

目前，关于中枢节点检测的方法都是针对静态功能脑网络设计的，很难捕获中枢节点随时间产生的动态变化，导致检测的结果不具有时间一致性且无法保证中枢节点的变化与功能脑网络中呈现的认知变化相一致。

发明内容

本发明提出一种基于图论的功能脑网络中枢节点的动态检测方法。针对上述静态检测方法的不足，本文在多变量中枢节点检测方法的基础上进行了改进，使得能够检测出更加可靠的符合神经科学认知活动的中枢节点。

首先，利用滑动窗口的技术，把血氧信号以时间为维度平均分成几段。在每一段时间的滑动窗口内，检测出对应时间窗口内的中枢节点，从而得到一条中枢节点随着滑动窗口移动的一条变化轨迹(如图1所示)。最后，把这条变化轨迹作为一种约束作用在多变量检测的方法上，从而能够动态的检测出更加可靠准确的中枢节点。

本发明解决其技术问题所采用的方案如下步骤：

步骤1:使用皮尔逊相关性在每个滑动窗口中构造功能脑网络；

基本数据格式和预备知识，用G＝(V,W)来表示一个功能脑网络，其中V表示N个脑区节点的集合，是一个N×N的邻接矩阵。w_ij表示邻接矩阵W中第i行第j列的元素，具体表示第i个脑区节点和第j个脑区节点的关联度；接下来计算拉普拉斯矩阵L＝D-W，其中D是一个对角元素为/>的度矩阵。

在图论中，图G由一组正交向量Φ形成，这些正交向量通过对拉普拉斯矩阵L进行特征分解来获得，即L＝Φ^TΛΦ，对角矩阵Λ＝diag[λ₁,λ₂,…,λ_N]，且λ₁,λ₂,…,λ_N表示以升序排序的特征值。

所述的特征值与正交向量一一对应，是拉普拉斯矩阵L的特征值；

如果所有脑区节点都连接在一个功能脑网络中，没有分离的部分，那么其中最小的特征值λ₁为零。零特征值的数量等于子图的数量。即当λ₁和λ₂都的值都为零的时候，子图的数量为2。

步骤2：改进的多变量中枢节点检测

2-1.功能脑网络中的每个脑区节点v_i都与一个选择向量中的二值标记s_i相关联，其中s_i＝0表示脑区节点v_i为中枢节点，s_i＝1表示不是中枢节点。通过移除选中的脑区节点后对功能脑网络的破坏程度，来判断是否是希望得到的中枢节点，其中破坏程度通过移除该脑区节点后剩余功能脑网络的零特征值的数量情况来判断。

判断准则如下：若剩余功能脑网络的零特征值的数量在删除选中的节点后，相对于删除别的节点增加的更多，则将该节点作为候选的中枢节点。

2-2.由于一次性需要选择多个中枢节点，所以对选择向量s＝[s₁,s₂,…,s_n]的估计是一个NP难的问题。根据Ky Fan定理，把问题转换成最小化前K个最小特征值的总和，

设功能脑网络的N个脑区节点中总共有K个中枢节点，则计算如下：

然后进一步推导得到最终的目标函数：

其中，λ_i表示第i个特征值；表示一个矩阵；

2-3.将选择向量s的每个元素放到对角矩阵S的对角线上。其中L_s＝D-S^TWS表示剩余功能脑网络的拉普拉斯矩阵，下标s是用于区分的标记；表示剩余功能脑网络中每个脑区节点的K维矩阵，并受F^TF＝I正交约束，因此需要优化求解F和S：

优化F：固定对角矩阵S，得到F的闭式解，闭式解是L_s的前K项最小特征值所对应的K个正交向量。

优化S：固定F，优化对角矩阵S的目标函数变为：

min_s trace(S^TAS) (3)

其中是一个元素a_ij＝w_ij||f_i-f_j||²的N×N的矩阵，且f_i和f_j分别表示F矩阵的第i行向量和第j行向量。

2-4.由于优化后的目标函数不是严格凸的，s是二值的选择向量不利于优化，因此引入辅助向量并将上述优化后的目标函数化简为:

其中P是由辅助向量p导出的另一个对角矩阵。直观地说，选择向量s是辅助向量p的二值化结果。

同样在方程(4)中交替求解P和S：

优化F：固定对角矩阵S，得到F的闭式解，闭式解是L_s的前K项最小特征值所对应的K个正交向量。

优化S和P：固定F，利用增广拉格朗日乘子法优化求解对角矩阵S和辅助矩阵P。

步骤3：功能脑网络中枢节点的动态检测

本发明的中枢节点检测对动态功能网络的算法框架如图2所示。

具体来如下所示：

3-1.将整个血氧信号分段成T组重叠的滑动窗口。在每个时刻t，估计每个脑区节点v_i对应的p_i，p_i是辅助向量p中的第i个元素。

3-2.连接p_i形成一条轨迹其中/>是连续函数p_i(τ)在时间点{τ_t}的离散采样。因此连续函数p_i(τ)在任意其它时间点τ的值都可以使用径向基函数(RBFs)计算：

其中参数σ用于控制轨迹平滑的强度，表示第i个脑区节点在时刻t的辅助向量的权值，即/>表示p_i在时刻t的权值；τ-τ_t表示时间差。

3-3.给定轨迹径向基函数集合/>通过以下方法计算：

其中，参数λ用于控制连续函数p_i(τ)的时间相关性的强度。

3-4.使用改进后由辅助向量构成的对角矩阵P，来初始化每个滑动窗口的选择向量s组成的对角矩阵S，如图2中的箭头所示。

3-5.将改进后由辅助向量构成的对角矩阵P带入步骤2-3重复优化求解直到目标收敛。

本发明有意效果如下：

本发明方法的输出是随着脑网络的变化而变化的中枢节点，这些中枢节点不仅在每个时间点上都比传统方法更精确，而且还保留了在时间点之间的相关性信息。

附图说明

图1是中枢节点动态检测方法的整体框架；

图2是中枢节点动态检测方法的算法框架；

图3是传统方法和本方法稳定性的分析图；

图4是传统方法和本方法对一个正常人中枢节点检测的可视化结果图；

图5是传统方法和本方法对一个强迫症患者中枢节点检测的可视化结果图；

具体实施方式

下面结合附图及实验，对本发明进行进一步的详细说明。

本发明利用滑动窗口的技术，把血氧信号以时间为维度平均分成几段。在每一段时间的滑动窗口内，检测出对应时间窗口内的中枢节点，从而得到一条中枢节点随着滑动窗口移动的一条变化轨迹(如图1所示)。最后，把这条变化轨迹作为一种约束作用在多变量检测的方法上，从而能够动态的检测出更加可靠准确的中枢节点。

本发明实施例实现的主要步骤：

步骤(1):选择63名正常人和62名强迫症患者的数据进行实验。每个受试者都有T1加权磁共振图像(具体参数为TR＝8毫秒，TE＝1.7ms毫秒，翻转角度＝20°，分辨率＝1.0×1.0×1.0mm²)和静息态功能磁共振数据(具体参数为TR＝2s，TE＝60ms毫秒，翻转角度＝90°，分辨率＝3.0×3.0×4.0mm²)，每个受试者都包含230个检测时间点。

步骤(2):将所有这些实验数据配准到AAL模板中，划分成116个脑区。并且计算这些脑区之间的相关性得到相应的功能脑网络连接矩阵W，作为实验的输入。

步骤(3):设置实验参数。

滑动窗口的大小设置成整个时间点的10％大小；中枢节点的检测数目设置为12；径向基RBF的参数σ设置成0.7；通过网格搜索法找到目标函数(2)的最优的λ参数为0.6。

步骤(4)：最后对目标函数依次优化F和S，直到收敛，得到最终的结果。

实验结果分析：

对于每一个个体，统计出每个滑动窗口内被检测为中枢节点的个数，构建相应的直方图，最后计算对应的熵值。

熵越低表明，中枢节点在整个检测过程中的变化也是少的，如图3所示，我们可以看到熵值都处于对角线以上，表明我们的方法检测出来的熵值更低，更加稳定。

我们进一步可视化了其中一段时间中枢节点的变化结果，如图4和图5所示,我们的方法检测出来的中枢节点在3个TR内(约6秒)没有发生变化，而传统方法检测出来的中枢节点在第二个TR时刻都发生了跳变，并在下一时刻恢复到之前的状态，这种情况是不寻常，不符认知变化的规律。因此，可以看出我们的方法比传统方法更加稳定和可靠。

Claims (4)

1.一种基于图论的功能脑网络中枢节点的动态检测方法，其特征在于利用滑动窗口的技术，把血氧信号以时间为维度平均分成几段；在每一段时间的滑动窗口内，检测出对应时间窗口内的中枢节点，从而得到一条中枢节点随着滑动窗口移动的一条变化轨迹；最后将这条变化轨迹作为一种约束作用在多变量检测的方法上，从而能够动态的检测出更加可靠准确的中枢节点；

包括如下步骤：

步骤1:使用皮尔逊相关性在每个滑动窗口中构造功能脑网络；

用G＝(V,W)来表示一个功能脑网络，其中V表示N个脑区节点的集合，是一个N×N的邻接矩阵；w_ij表示邻接矩阵W中第i行第j列的元素，具体表示第i个脑区节点和第j个脑区节点的关联度；接下来计算拉普拉斯矩阵L＝D-W，其中D是一个对角元素为的对角矩阵；

在图论中，图G由一组正交向量Φ形成，这些正交向量通过对拉普拉斯矩阵L进行特征分解来获得，即L＝Φ^TΛΦ，对角矩阵Λ＝diaag[λ₁,λ₂,…,λ_N]，且λ₁,λ₂,…,λ_N表示以升序排序的特征值；

步骤2：改进的多变量中枢节点检测

2-1.功能脑网络中的每个脑区节点v_i都与一个选择向量中的二值标记s_i相关联，其中s_i＝0表示脑区节点v_i为中枢节点，s_i＝1表示不是中枢节点；通过移除选中的脑区节点后对功能脑网络的破坏程度，来判断是否是希望得到的中枢节点，其中破坏程度通过移除该脑区节点后剩余功能脑网络的零特征值的数量情况来判断；

2-2.由于一次性需要选择多个中枢节点，所以对选择向量s＝[s₁,s₂,…,s_n]的估计是一个NP难的问题；根据Ky Fan定理，把问题转换成最小化前K个最小特征值的总和；

设功能脑网络的N个脑区节点中总共有K个中枢节点，则计算如下：

然后进一步推导得到最终的目标函数：

其中，λ_i表示第i个特征值；表示剩余功能脑网络中每个脑区节点的K维矩阵；

2-3.将选择向量s的每个元素放到对角矩阵S的对角线上；其中L_s'＝D-S^TWS表示剩余功能脑网络的拉普拉斯矩阵，下标s'是用于区分的标记；表示剩余功能脑网络中每个脑区节点的K维矩阵，并受F^TF＝I正交约束，因此需要优化求解F和S：

优化F：固定对角矩阵S，得到F的闭式解，闭式解是L_s'的前K项最小特征值所对应的K个正交向量；

优化S：固定F，优化对角矩阵S的目标函数变为：

其中是一个元素a_ij＝w_ij||f_i-f_j||²的N×N的矩阵，且f_i和f_j分别表示F矩阵的第i行向量和第j行向量；

2-4.由于优化后的目标函数不是严格凸的，s是二值的选择向量不利于优化，因此引入辅助向量并将上述优化后的目标函数化简为:

其中P是由辅助向量p导出的另一个对角矩阵即辅助矩阵；直观地说，选择向量s是辅助向量p的二值化结果；

同样在方程(4)中交替求解P和S：

优化F：固定对角矩阵S，得到F的闭式解，闭式解是L_s'的前K项最小特征值所对应的K个正交向量；

优化S和P：固定F，利用增广拉格朗日乘子法优化求解对角矩阵S和辅助矩阵P；

步骤3：功能脑网络中枢节点的动态检测

3-1.将整个血氧信号分段成T组重叠的滑动窗口；在每个时刻t，估计每个脑区节点v_i对应的p_i，p_i是辅助向量p中的第i个元素；

3-2.连接p_i形成一条轨迹其中/>是连续函数p_i(τ)在时间点{τ_t}的离散采样；因此连续函数p_i(τ)在任意其它时间点τ的值都可以使用径向基函数计算：

其中参数σ用于控制轨迹平滑的强度，表示第i个脑区节点在时刻t的辅助向量的权值，即/>表示p_i在时刻t的权值；τ-τ_t表示时间差；

3-3.给定轨迹径向基函数集合/>通过以下方法计算：

其中，参数λ用于控制连续函数p_i(τ)的时间相关性的强度；

3-4.使用辅助矩阵P，来初始化每个滑动窗口的选择向量s组成的对角矩阵S；

3-5.将辅助矩阵P带入步骤2-3重复优化求解直到目标收敛。

2.根据权利要求1所述的一种基于图论的功能脑网络中枢节点的动态检测方法，其特征在于所述的特征值与正交向量一一对应，是拉普拉斯矩阵L的特征值。

3.根据权利要求1或2所述的一种基于图论的功能脑网络中枢节点的动态检测方法，其特征在于如果所有脑区节点都连接在一个功能脑网络中，没有分离的部分，那么其中最小的特征值λ₁为零；零特征值的数量等于子图的数量；即当λ₁和λ₂都的值都为零的时候，子图的数量为2。

4.根据权利要求3所述的一种基于图论的功能脑网络中枢节点的动态检测方法，其特征在于2-1中判断的准则如下：若剩余功能脑网络的零特征值的数量在删除选中的节点后，相对于删除别的节点增加的更多，则将该选中的节点作为候选的中枢节点。