基于网络特征熵定量刻画脑功能网络的方法

摘要

本发明公开了一种基于网络特征熵定量刻画脑功能网络的方法，步骤包括：采集静息状态下的功能磁共振成像数据并进行预处理，选择包含大脑皮层功能分区的脑模板，将脑模板中的单个脑区看成脑功能网络中的一个节点，提取功能磁共振成像数据在每个节点的平均时间序列信号，计算节点之间的功能连接，得到功能连接矩阵，将所述功能连接矩阵进行二值化得到邻接矩阵，形成脑功能网络，计算脑功能网络中每一个节点的特征向量中心度及能量集中度，根据脑功能网络中每一个节点的能量集中度计算大脑网络特征熵。本发明通过大脑网络特征熵定量刻画脑功能网络以揭示大脑发育及老化规律，具有原理简单、实现方便、被试内结果稳定、不同类别样本灵敏性好的优点。

说明书

技术领域

本发明涉及基于功能磁共振成像(fMRI)数据的大脑网络研究领域，具体涉及一种基于网络特征熵定量刻画脑功能网络的方法。

背景技术

功能磁共振成像的原理是通过磁振造影量测神经元活动引起的血液动力的改变。例如：刺激某些感官时，大脑皮层具体的某些脑功能区激活，这些激活的脑区发出神经信号，利用磁振造影获得神经信号，得到一系列图像。人们执行运动、感觉、高级认知等各种复杂的任务时，会激起脑内某些神经信号，这种生理性的大脑活动引起局部脑血流、脑血容、能量代谢的变化，从而导致组织磁性的变化。fMRI作为一种活体脑功能成像技术，因其具有非侵入性，理想时空分辨率等优点，迅速成为脑神经科学研究最常用的脑成像技术。

不管是生物学领域、技术领域、社会领域或者其他的科学领域都产生了类似地相互有联系的数据集。在过去的几十年，为了刻画这些数据集，复杂网络分析应运而生。复杂网络分析主要通过图论定量刻画给定网络的拓扑特性。大脑连接数据是指有功能连接或结构追踪连接形成的网络，它与生物或者物理系统有共同的特征，同时研究者发现大脑从微小的神经元到脑区等不同的层次上都可以建立网络模型，这样利用复杂网络分析方法研究大脑网络属性为人们从整体上认识和研究大脑功能活动提供了一种全新的研究思路。

2006年，Achard Sophi等人利用离散小波变化转化fMRI数据估计刻画功能连接的频率相关性矩阵，得到结论：fMRI数据中低频振荡数据集具有小世界(small world)性质。2008年van den Heuvel MP等人根据静息功能核磁共振(rs-fMRI)数据基于体素(voxel-wise)方法，构建体素之间的功能连接，然后计算网络图的统计属性发现脑功能网络的聚类系数比对应随机网络的要大，并且平均路径长度要比随机网络短，表明了脑功能网络的小世界特性；另外，体素与体素之间功能连接的分布服从幂率指数近似等于2的幂率分布，表明脑功能网络的无标度特性(scale free network)。但是，以上得到的这些拓扑特性只刻画了脑功能网络的一个特征，这些特征是孤立的，它的应用受到局限。比如我们探索受到很多因素影响的大脑发育，必须将这些特性进行综合分析，提取网络的结构趋势特征，如何通过综合分析定量刻画脑功能网络以探索网络特征熵与大脑发育、老化规律，已经成为一项亟待解决的关键技术问题。

发明内容

本发明要解决的技术问题：针对现有技术的上述问题，提供通过大脑网络特征熵揭示大脑发育及老化规律，具有原理简单、实现简便、被试内结果稳定性高、不同类别灵敏性好的基于网络特征熵定量刻画脑功能网络的方法。

为了解决上述技术问题，本发明采用的技术方案为：

一种基于网络特征熵定量刻画脑功能网络的方法，步骤包括：

1)采集静息状态下的功能磁共振成像数据；

2)针对采集到的功能磁共振成像数据进行预处理；

3)选择包含大脑皮层功能分区的脑模板；

4)针对预处理后的功能磁共振成像数据，将脑模板上的不同脑区看成脑功能网络中的节点，从而提取每个节点的平均时间序列信号，计算节点之间的相关系数，得到功能连接矩阵，将所述功能连接矩阵进行二值化得到邻接矩阵，从而形成脑功能网络，所述邻接矩阵形成的脑功能网络刻画了节点间相互作用关系，反映了脑功能网络的拓扑属性；

5)计算脑功能网络中每一个节点的特征向量中心度及能量集中度；

6)根据脑功能网络中每一个节点的能量集中度计算大脑网络特征熵。

优选地，所述步骤4)的详细步骤包括：

4.1)针对预处理后的功能磁共振成像数据，将脑模板上的不同脑区看成脑功能网络中的节点，提取每个节点的平均时间序列信号；

4.2)采用式(1)计算节点之间的相关系数，得到一个对称的大小为n*n的时间序列相关矩阵A(a_ij)，其中n表示脑功能网络节点个数，a_ij表示节点i的平均时间序列信号x_i与节点j的平均时间序列信号x_j之间的相关系数，以相关系数表示节点之间的功能连接，因此所述时间序列相关矩阵A(a_ij)表示了脑功能网络中各节点间功能连接，称为功能连接矩阵；

$a_{ij}=\frac{E\left(x_i{x}_j\right)-E\left(x_i\right)E\left(x_j\right)}{\sqrt{E\left(x_i^2\right)-E^2\left(x_i\right)}\sqrt{E\left(x_j^2\right)-E^2\left(x_j\right)}}---\left(1\right)$

式(1)中，E(x_ix_j)表示x_ix_j数学期望，E(x_i)表示x_i的数学期望，E(x_j)表示x_j的数学期望，表示x_i²的数学期望，E²(x_i)表示x_i数学期望的平方，表示x_j²的数学期望，E²(x_j)表示x_j数学期望的平方；

4.3)针对所述功能连接矩阵，首先根据任意两节点之间功能连接对应的统计值p，通过Bonferroni校正，将对应统计值p大于预设统计值的节点对之间的功能连接设为0，将对应统计值p小于或等于预设统计值的节点对之间的功能连接设为1，其次将所述功能连接矩阵中小于0的元素设为0，从而将所述功能连接矩阵二值化得到邻接矩阵形成脑功能网络，所述脑功能网络刻画了节点间相互作用关系，反映了脑功能网络的拓扑属性。

优选地，所述步骤4.3)中的预设统计值为0.05。

优选地，所述步骤5)的详细步骤包括：

5.1)根据式(2)计算所述邻接矩阵的最大特征值对应特征向量的第i个分量作为邻接矩阵中第i个节点的特征向量中心度；

$e_i=\frac{1}{\lambda }\underset{j=1,j\ne i}{\overset{n}{\Sigma }}{a}_{ij}{e}_j---\left(2\right)$

式(2)中，e_i表示邻接矩阵中第i个节点的特征向量中心度，a_ij表示邻接矩阵中第i个节点与第j个节点之间的功能连接，λ表示邻接矩阵对应的最大特征值，e_j表示λ对应的特征向量，n表示脑功能网络节点个数；

5.2)根据式(3)计算邻接矩阵中每一个节点的能量集中度；

$I_i=\frac{e_i}{\underset{i=1}{\overset{n}{\Sigma }}{e}_i}---\left(3\right)$

式(3)中，I_i表示邻接矩阵中第i个节点的能量集中度，e_i表示邻接矩阵中第i个节点的特征向量中心度，n表示脑功能网络节点个数。

优选地，所述步骤6)中计算大脑网络特征熵的函数表达式如式(4)所示；

$E=-\underset{i=1}{\overset{n}{\Sigma }}{I}_i{\mathrm{lnI}}_i---\left(4\right)$

式(4)中，E表示大脑网络特征熵，I_i表示邻接矩阵中第i个节点的能量集中度，n表示脑功能网络节点个数。

优选地，所述步骤2)的详细步骤包括：

2.1)去除功能磁共振成像数据的前指定数量帧的图像；

2.2)剔除功能磁共振成像数据中在头部挪动幅度超过指定阈值的被试；

2.3)针对功能磁共振成像数据进行层间时间校正；

2.4)针对功能磁共振成像数据进行头动校正；

2.5)针对功能磁共振成像数据进行空间标准化处理；

2.6)利用指定大小半高宽的高斯核函数对功能磁共振成像数据进行空间平滑处理；

2.7)利用指定频率区间对功能磁共振成像数据进行时域滤波；

2.8)针对功能磁共振成像数据进行回归全脑平均信号、白质、脑积液平均信号和头动参数。

优选地，所述步骤2.6)中指定大小半高宽具体是指6mm的半宽高。

优选地，所述步骤2.7)中指定频率区间具体是指0.01～0.08Hz。

本发明基于网络特征熵定量刻画脑功能网络的方法具有下述优点：

1、本发明将脑模板中的单个脑区看成网络中的节点，通过提取预处理后功能磁共振成像数据在每个节点的平均时间序列信号，计算节点之间的功能连接，得到功能连接矩阵，将功能连接矩阵进行二值化得到邻接矩阵，形成脑功能网络，针对脑功能网络的邻接矩阵，分别计算每一个节点的特征向量中心度及能量集中度，挖掘出了与网络本质拓扑属性相关的特征。

2、本发明通过已有的香农信息熵理论，根据每一个节点的能量集中度计算大脑网络特征熵，由此可以通过比较计算不同发育阶段的大脑网络，探索大脑发育及老化规律，通过重复实验验证了网络特征熵在个体上的稳定性，为利用网络特征熵探索大脑发育、老化规律提供方法和思路，具有原理简单、实现简便、被试内结果稳定性高、不同类别灵敏性好的优点。

附图说明

图1为本发明实施例方法的基本流程示意图。

图2为本发明实施例中被试年龄分布图。

具体实施方式

如图1所示，本实施例基于网络特征熵定量刻画脑功能网络的方法的步骤包括：

1)采集静息状态下的功能磁共振成像数据；

2)针对采集到的功能磁共振成像数据进行预处理；

3)选择包含大脑皮层功能分区的脑模板；

4)针对预处理后的功能磁共振成像数据，将脑模板上的不同脑区看成脑功能网络中的节点，从而提取每个节点的平均时间序列信号，计算节点之间的相关系数，得到功能连接矩阵，将所述功能连接矩阵进行二值化得到邻接矩阵，从而形成脑功能网络，所述邻接矩阵形成的脑功能网络刻画了节点间相互作用关系并反映了脑功能网络的拓扑属性；

5)计算脑功能网络中每一个节点的特征向量中心度及能量集中度；

6)根据脑功能网络中每一个节点的能量集中度计算大脑网络特征熵。

本实施例步骤1)中，采集静息状态(resting state)下的功能磁共振成像(functional magnetic resonance imaging，简称fMRI)数据时，使用的设备为核磁共振扫描仪，且选择快速梯度回波方式，扫描前需要设定扫描周期、层厚、层数、扫描矩阵、扫描视野等扫描参数，扫描过程中要求被试保持头部不动，闭上眼睛放松，但保持清醒。

本实施例中，步骤2)的详细步骤包括：

2.1)去除功能磁共振成像数据的前指定数量帧的图像，用于以消除磁场饱和给数据带来的不准确影响，本实施例中具体为去掉前5帧图像；

2.2)剔除功能磁共振成像数据中在头部挪动幅度超过指定阈值的被试；

2.3)针对功能磁共振成像数据进行层间时间校正；

2.4)针对功能磁共振成像数据进行头动校正；

2.5)针对功能磁共振成像数据进行空间标准化处理；

2.6)利用指定大小半高宽的高斯核函数对功能磁共振成像数据进行空间平滑处理；本实施例中，步骤2.6)中指定大小半高宽具体是指6mm的半宽高；

2.7)利用指定频率区间对功能磁共振成像数据进行时域滤波；本实施例中，步骤2.6)中指定频率区间具体是指0.01～0.08Hz；

2.8)针对功能磁共振成像数据进行回归全脑平均信号、白质、脑积液平均信号和头动参数，以降低硬件噪声。

本实施例中，步骤2.2)～2.8)具体是使用Statistical Parametric Mapping软件(SPM8)在Matlab运行下实现，通过上述处理，能够提高功能磁共振成像数据的准确度。

脑模板是前人通过任务态激活区检测等方法确定的大脑皮层功能分区，本实施例步骤3)中，选择包含大脑皮层功能分区的脑模板具体为Montreal Neurological Institute(MNI)机构提供的Anatomical Automatic Labeling(AAL)模板，它一共包含116个脑区；此外，也可以根据需要选择其他公认较稳定的脑模板。

本实施例中，步骤4)的详细步骤包括：

所述步骤4)的详细步骤包括：

4.1)针对预处理后的功能磁共振成像数据，将脑模板上的不同脑区看成脑功能网络中的节点，提取每个节点的平均时间序列信号；

4.2)采用式(1)计算节点之间的相关系数，得到一个对称的大小为n*n的时间序列相关矩阵A(a_ij)，其中n表示脑功能网络节点个数，a_ij表示节点i的平均时间序列信号x_i与节点j的平均时间序列信号x_j之间的相关系数，以相关系数表示节点之间的功能连接，因此所述时间序列相关矩阵A(a_ij)表示了脑功能网络中各节点间功能连接，称为功能连接矩阵；

$a_{ij}=\frac{E\left(x_i{x}_j\right)-E\left(x_i\right)E\left(x_j\right)}{\sqrt{E\left(x_i^2\right)-E^2\left(x_i\right)}\sqrt{E\left(x_j^2\right)-E^2\left(x_j\right)}}---\left(1\right)$

式(1)中，E(x_ix_j)表示x_ix_j数学期望，E(x_i)表示x_i的数学期望，E(x_j)表示x_j的数学期望，表示x_i²的数学期望，E²(x_i)表示x_i数学期望的平方，表示x_j²的数学期望，E²(x_j)表示x_j数学期望的平方；

4.3)针对所述功能连接矩阵，首先根据任意两节点之间功能连接对应的统计值p，通过Bonferroni校正，将对应统计值p大于预设统计值的节点对之间的功能连接设为0，将对应统计值p小于或等于预设统计值的节点对之间的功能连接设为1，其次将所述功能连接矩阵中小于0的元素设为0，从而将所述功能连接矩阵二值化得到邻接矩阵形成脑功能网络，所述脑功能网络刻画了节点间相互作用关系，反映了脑功能网络的拓扑属性。

本实施例中，步骤4.3)中的预设统计值为0.05。

本实施例中，步骤5)的详细步骤包括：

5.1)根据式(2)计算所述邻接矩阵的最大特征值对应特征向量的第i个分量作为邻接矩阵中第i个节点的特征向量中心度(eigenvector centrality)；

$e_i=\frac{1}{\lambda }\underset{j=1,j\ne i}{\overset{n}{\Sigma }}{a}_{ij}{e}_j---\left(2\right)$

式(2)中，e_i表示邻接矩阵中第i个节点的特征向量中心度，a_ij表示邻接矩阵中第i个节点与第j个节点之间的功能连接，λ表示邻接矩阵对应的最大特征值，e_j表示λ对应的特征向量，n表示脑功能网络节点个数；

5.2)根据式(3)计算邻接矩阵中每一个节点的能量集中度；

$I_i=\frac{e_i}{\underset{i=1}{\overset{n}{\Sigma }}{e}_i}---\left(3\right)$

式(3)中，I_i表示邻接矩阵中第i个节点的能量集中度，e_i表示邻接矩阵中第i个节点的特征向量中心度，n表示脑功能网络节点个数；

本实施例中，步骤6)中计算大脑网络特征熵的函数表达式如式(4)所示；

$E=-\underset{i=1}{\overset{n}{\Sigma }}{I}_i{\mathrm{lnI}}_i---\left(4\right)$

式(4)中，E表示大脑网络特征熵，I_i表示邻接矩阵中第i个节点的能量集中度，n表示脑功能网络节点个数。式(4)所示的函数表达式为通过已有的香农的信息熵理论大脑网络的熵，基于网络特征熵定量刻画脑功能网络，可以用于探索大脑发育、老化规律。

为了评估大脑网络特征熵在脑功能网络应用的性能，二个因素必须被考虑：1、重测信度(test-retest reliability)：重测信度的高低反映了被试内网络特征熵可变性大小，衡量了网络特征熵在时域上的稳定性，2、灵敏性：通过双样本t检验或非参数检验检验不同类别(例如：男女之间)被试之间的差异是否显著。本实施例针对性别匹配共97人(男性54，女性43)的同一批被试，通过得到不同时间情况下采集到的两个Scan fMRI数据(Scan之间相隔一周)，计算网络特征熵的重测信度和灵敏性。被试平均年龄为27.83岁，方差为5.09岁，年龄分布如图2所示。

重测信度是指同一个静息态功能磁共振实验在不同的时间对同一被试群体施测所得到的测量结果的一致性程度。在功能连接组学中，常用组内相关系数intra-class correlation(ICC)定量估计一个度量的重测信度。参考Zuo Xi nian 2014年发表综述文章test-retest reliabilities of resting-state FMRI measurements inhuman brain functional connectomics:A systems neuroscience perspective中将重测信度分为五个等级：较低：0<ICC≤0.2，一般：0.2<ICC≤0.4，中等：0.4<ICC≤0.6，显著：0.6<ICC≤0.8，较高：0.8<ICC≤1。本实施例中，利用组内相关系数计算大脑网络特征熵的重测信度的步骤包括：

A1)将两个Scan fMRI数据(Scan1和Scan2)分别采用本实施例基于网络特征熵定量刻画脑功能网络的方法计算大脑网络特征熵，被试个数n＝97，测量次数d＝2，网络特征熵为E，E_ij(i＝1,…d；j＝1,…,n)是指第j个被试在第i次测量时得到的网络特征熵。将计算得到的大脑网络特征熵利用线性混合模型Linear>ij＝λ_0j+e_ij,λ_0j＝μ₀₀+p_0j，其中μ₀₀是一个固定的参数，它是指所有被试d次测量的网络特征熵的均值，λ_0j为中间参数，p_0j和e_ij分别是服从均值为0方差为和的独立正态分布随机误差，p_0j指的是被试间的误差，e_ij指的是测量误差。如果μ₀₀显著不为0，则大脑网络特征熵的组内平均系数ICC(E)定义如式(5)所示。

$ICC\left(E\right)=\frac{{\mathrm{MS}}_b\left(E\right)-{\mathrm{MS}}_w\left(E\right)}{{\mathrm{MS}}_b\left(E\right)+(2-1){\mathrm{MS}}_w\left(E\right)}=\frac{{\sigma }_p^2}{{\sigma }_p^2+{\sigma }_e^2}---\left(5\right)$

式(5)中，MS_b(E)为两两样本间方差的总和，MS_w(E)为所有被试内均方误差，和为方差。

A2)计算SS_p和SS_e。

根据线性混合模型E_ij，为了估计组内相关系数ICC，根据式(6)和(7)计算SS_p和SS_e。

${\mathrm{SS}}_p=2\underset{i=1}{\overset{n}{\Sigma }}{(\overline{Y_{i.}}-\overline{Y_{\mathrm{..}}})}^2---\left(6\right)$

${\mathrm{SS}}_e=\underset{i=1}{\overset{n}{\Sigma }}\underset{j=1}{\overset{2}{\Sigma }}{(\overline{Y_{ij}}-\overline{Y_{i.}}-\overline{Y_{.j}}+\overline{Y_{\mathrm{..}}})}^2---\left(7\right)$

式(6)和(7)中，SS_p为组间差平方和，SS_e为总误差平方和，为第j个样本在第i次测量时得到的网络特征熵，n为样本个数，其中第i个样本组间样本均值、第j组组内样本均值、总样本均值三者的表达式详见式(8)。

$\left(\begin{array}{c}\overline{Y_{i.}}=\frac{\underset{j=1}{\overset{n}{\Sigma }}{Y}_{ij}}{n}\\ \overline{Y_{.j}}=\frac{\underset{i=1}{\overset{d}{\Sigma }}{Y}_{ij}}{d}\\ \overline{Y_{\mathrm{..}}}=\frac{\underset{i=1}{\overset{d}{\Sigma }}\underset{j=1}{\overset{n}{\Sigma }}{Y}_{ij}}{nd}\end{array}\right)---\left(8\right)$

式(8)中，d为测量次数，其余参数含义与式(6)和(7)相同，在此不再赘述。

A3)估计ICC(E)。

MS_b＝SS_p/(n-1)和MS_w＝SS_e/(n-1)(d-1)的数学期望分别是

因此，和的无偏估计如式(9)所示；

$\left(\begin{array}{c}{\hat{\sigma }}_p^2=\frac{{\mathrm{MS}}_b-{\mathrm{MS}}_w}{2}\\ {\hat{\sigma }}_e^2={\mathrm{MS}}_w\end{array}\right)---\left(9\right)$

式(9)中，为方差的无偏估计，为方差的无偏估计。

那么ICC(E)的估计如式(10)所示；

$ICC\left(E\right)=\frac{{\hat{\sigma }}_p^2}{{\hat{\sigma }}_p^2+{\hat{\sigma }}_e^2}=\frac{{\mathrm{MS}}_b-{\mathrm{MS}}_w}{{\mathrm{MS}}_b+(d-1){\mathrm{MS}}_w}---\left(10\right)$

式(10)中，为方差的无偏估计，为方差的无偏估计。本实施例中，最终计算得到大脑网络特征熵的重测信度ICC(E)的估计值为0.96，本实施例基于网络特征熵定量刻画脑功能网络的方法的重测信度值属于较高等级，本实施例基于网络特征熵定量刻画脑功能网络的方法中的大脑网络特征熵这个度量被试内可变性小，时域稳定性高。

灵敏性是指同一个静息态功能磁共振实验对不同类别的被试群体施测所得到的测量结果的差异性程度。在统计学中，我们常采用双样本t检验检验两组不同群体的差异性，显著性p值设置为0.05，若p<0.05，则两组不同群体在该度量下具有显著性差异，否则两组不同群体在该度量下没有显著差异。本实施例中，确定灵敏性的步骤包括：

B1)同一个被试scan1和scan2时间序列合在一起，计算每个被试的网络特征熵。

B2)建立线性回归模型，E＝β₀+β₁×age，其中E表示可观测的变量具体指的是网络特征熵，得到的残差用于后续分析，其中β₀、β₁表示回归系数，age为自变量。

B3)将得到的残差根据性别分成两组，

B4)利用双样本t检验，检验男女两组被试网络特征熵的差异，计算得到p＝0.024，根据显著性水平p＝0.05知p<0.05。

综上所述，目前对于大脑网络的研究多为定性分析，本实施例基于网络特征熵定量刻画脑功能网络的方法结合信息熵的定义，探究网络特征熵的定义，给出能刻画大脑网络复杂程度的熵定义，并从考察度量个体内的稳定性(重测信度)和区分类别的灵敏性验证了网络特征熵的可靠性，进而比较发育、老化不同时期大脑的网络特征熵，探寻网络特征熵与大脑发育、老化规律之间存在的关系，以及发育、老化速率，从而研究大脑发育、老化规律。

以上所述仅是本发明的优选实施方式，本发明的保护范围并不仅局限于上述实施例，凡属于本发明思路下的技术方案均属于本发明的保护范围。应当指出，对于本技术领域的普通技术人员来说，在不脱离本发明原理前提下的若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。