CN110889501B - 一种基于非平稳动态贝叶斯网的人脑效应连接网络构建方法 - Google Patents

Abstract

本发明公开了一种基于非平稳动态贝叶斯网的动态人脑效应连接识别方法，获取功能磁共振成像数据；使用DPABI软件包对原始图像数据进行数据预处理；选取需要识别效应连接的脑区作为感兴趣区域；使用非平稳动态贝叶斯网方法，对提取的感兴趣区域进行脑效应连接的识别；观察由非平稳动态贝叶斯网构建的不同时刻下的动态脑效应连接网络，对不同时刻的脑效应连接网络进行分析，从而进一步理解大脑的内在机理和运行模式。本发明将动态贝叶斯网方法作为基础框架方法，该方法可以从fMRI时间序列数据中，构建出随时间变化的动态人脑效应连接网络。

Description

一种基于非平稳动态贝叶斯网的人脑效应连接网络构建方法

技术领域

本发明涉及功能磁共振成像数据的动态人脑效应连接网络构建方法，特别是一种基于非平稳动态贝叶斯网的方法。

背景技术

人脑是目前所知宇宙间最为复杂的系统之一，为了系统全面准确地探索人脑的工作机制，脑科学的研究须从连接和网络的层面进行。功能磁共振成像技术(functionalmagenetic resonance imaging,fMRI)是一种非侵入性活体脑功能成像技术。由于它具有可靠的理论基础、较高的时空分辨率，为认知神经科学的实验研究提供了有利条件，同时又为理解人类大脑提供了强有力的手段。特别是从fMRI数据中构建人脑效应连接网络，可以帮助了解复杂的人类大脑的运作，并为理解神经退行性疾病，如阿尔茨海默病(Alzheimersdisease，AD)、帕金森病和癫痫等提供帮助。具体来说，脑效应连接网络刻画了一个脑区施加于另一个脑区神经活动的因果效应，可以用一个有向图表示。其中，有向弧表示效应连接，弧的参数表示效应连接强度。

近年来，基于概率论、统计学和图论的贝叶斯网络(Bayesian network，BN)方法逐渐成为构建脑效应连接的一种重要而实用的方法。然而，贝叶斯网方法构建的人脑效应连接网络在时间上是静态的，且无法从fMRI时间序列数据中获取数据的时间信息。此外，贝叶斯网络是是一种有向无环图(Directed Acyclic Graph，DAG)，因此无法对循环或双向效应连接进行建模。为了克服贝叶斯网的这些问题，动态贝叶斯网络(Dynamic Bayesiannetwork，DBN)方法被提出，该方法能够对脑区之间的动态关系进行建模，有效地从fMRI时间序列数据中获取时间信息。

然而，目前对贝叶斯网方法和动态贝叶斯网方法的研究很大程度上依赖于一个假设，即整个fMRI时间序列的连接模式是静态的或假设时间序列是平稳，然而这种假设并不能真正反映人脑效应连接网络的结构。在这种假设下，人脑效应连接网络由大脑各区域之间的空间和时间恒定关系来表示，并侧重于在整个时间序列中评估统计依赖性的度量。另一方面，脑区的之间的相互作用和信号转导过程通常是随着外界刺激而改变的。例如，大脑进行不同任务或着疲劳的时候，都可能引发大脑不同的响应。因此，有必要对动态的脑效应连接网络进行构建，以反映大脑实际的运作机理。

发明内容

针对上述人脑效应连接网络构建方法存在的不足，本发明将动态贝叶斯网方法作为基础框架，提出了一种利用非平稳动态贝叶斯网络构建动态人脑效应连接网络的新方法，该方法可以从fMRI时间序列数据中，构建出随时间变化的动态人脑效应连接网络。

实现本发明的主要思路是：获取功能磁共振成像数据；使用DPABI(a toolbox forData Processing&Analysis for Brain Imagin)软件包对原始图像数据进行数据预处理；选取需要识别效应连接的脑区作为感兴趣区域(Region of Interest，ROI)；使用非平稳动态贝叶斯网方法，对提取的感兴趣区域进行脑效应连接的识别；观察由非平稳动态贝叶斯网构建的不同时刻下的动态脑效应连接网络，对不同时刻的脑效应连接网络进行分析，从而进一步理解大脑的内在机理和运行模式。

一种基于非平稳动态贝叶斯网的动态人脑效应连接识别方法包括以下步骤：

步骤1数据获取：本方法所使用的数据，包含仿真fMRI数据和真实fMRI数据，仿真fMRI数据用于检测方法性能，真实fMRI数据用于探究实际应用价值。仿真fMRI数据由于已知标准网络结构，因此可以用于检测算法有效性，衡量算法的性能。真实数据用于检测算法在真实fMRI数据情况中的动态脑效应连接识别能力。仿真fMRI数据用开源工具包SimTB进行生成。真实fMRI数据使用ADNI数据库中的fMRI数据。

步骤2数据预处理：真实fMRI数据使用DPABI软件包来进行数据预处理。

步骤3选取感兴趣区域(ROI)：针对仿真fMRI数据，从DPABI软件组件模板中选择10个脑区作为感兴趣区域。对于真实fMRI数据采用最传统的AAL模板，然后提取42个与AD相关脑区的体素时间序列。

步骤4数据离散化。由于网络结构使用贝叶斯评分(BDe,Bayesian Dirichletequivalent score)进行衡量，因此需要对连续的fMRI时间序列进行离散化处理。

步骤5使用非平稳动态贝叶斯网方法进行动态人脑效应连接网络构建。具体包括以下几个步骤：初始化参数，包括状态个数s，转移时间集t，最大父母节点个数m和马尔科夫阶数(延迟)l；初始化原始时刻网络，使用马尔科夫链蒙特卡洛方法搜索初始网络中贝叶斯评分最高的贝叶斯网络；确定转移集，寻找下一时刻效应连接发生变化的节点；网络动态转换，将初始网络与转移集进行结合，确定下一时刻的贝叶斯网；重复进行网络转移过程，直到网络转移到达最后时刻；确定最终贝叶斯评分最高的动态贝叶斯网，包含一个初始人脑效应连接网络和(s-1)个后续人脑效应连接网络，这s个网络共同组成一个完整的、随时间变化的动态人脑效应连接网络。

步骤6结果分析：对于构建的动态人脑效应连接网络结构(包含s个动态贝叶斯网络)，节点表示脑区(ROI)，有向弧为脑区之间的效应连接。在仿真fMRI数据集上，将本方法与传统贝叶斯网方法和动态贝叶斯网方法构建的脑效应连接网络进行对比，从而验证方法有效性。在真实fMRI数据中，使用本方法，对在不同年龄时的fMRI数据进行脑效应连接网络构建，通过分析网络随时间变化的差异，寻找年龄增长而变化的因素。

与现有技术相比，本发明具有以下明显的优势和有益效果；

(1)本发明提出了一种非平稳动态贝叶斯网方法，可以将脑效应连接网络随时间的变化刻画出来。

(2)本发明使用非平稳动态贝叶斯网进行动态脑效应连接网络结构学习，相比传统的动态贝叶斯网的方法，新方法不仅能够捕获fMRI数据中的时间信息，而且能够获得网络结构随时间变化的动态人脑效应连接网络。

(3)本发明通过数据离散方法，去除方法对数据分布的依赖性，相较其他算法，在仿真数据集上获得了更准确的网络。

(4)本发明在真实fMRI数据集上学习的脑效应连接网络结构合理可靠，可以观察脑效应连接网络结构随年龄增大而变化的动态过程，具有实际的意义和价值，可为脑网络异常检测提供帮助。

附图说明

图1为本发明所涉及方法的流程图。

图2仿真fMRI数据生成示意图。

图3仿真fMRI数据感兴趣区域选择。

图4非平稳动态贝叶斯网构建过程。

图5实施例的AD被试动态脑效应连接网络结构连接图。

具体实施方式

下面阐述本发明的具体实施方式和详细步骤，本发明所涉及方法的流程图如图1所示，具体包括：

(步骤1)数据获取。

首先为了验证本发明方法的有效性，使用开源的仿真fMRI数据生成工具包(Asimulation toolbox for fMRI data，SimTB)生成一组基于事件相关的仿真数据。工具箱允许在时空可分离性模型下灵活地生成功能磁共振成像数据集，并被设计用于各种分析方法的测试。在SimTB中，可以完全控制数据生成，包括创建和操作空间源、实现与块和事件相关的实验设计、包含组织特定基线、模拟头部运动等。具体来说，首先确定被试者个数(20人)，然后确定扫描时间即时间点个数(180个)，接着在给定的脑区模板中选择相应的感兴趣区域(10个)。设置特定的任务事件，定义一个初始静息态和3个任务状态(共计3个状态)。对模板的体素数据进行平移，旋转，调整大小，规格化到空间图中。最后加入噪声，并与血液动力学相应函数做卷积操作。最终将得到10个脑区，180个时间点，3个状态的20名独立被试的体素时间序列数据。仿真数据生成示意图如图2所示。

真实fMRI数据来自于ADNI(Alzheimer's Disease Neuroimaging Initiative)数据库，ADNI数据库是一个公共数据库，由超过1500名被试的采集数据组成，被试的年龄在55到100岁之间，包括早期轻度认知障碍(Early Mild Cognitive Impairment，EMCI)被试，晚期轻度认知障碍(Late Mild Cognitive Impairment，LMCI)被试，阿兹海默症(Alzheimer's disease，AD)被试以及健康对照组(Health Control，HC)被试等。本发明拟观察AD被试随着年龄增长，脑效应连接网络的变化情况。故选取了一名AD被试从86岁到89岁采集的多个fMRI数据作为一个数据集。采集fMRI数据的具体参数设置如下：西门子3T超导型磁共振仪,标准8通道头线圈。重复时间为3000ms，回波时间为30ms，翻转角为90°，矩阵为64×64，层厚为4mm，共48层，被试采集的时间点为160个时间点。

步骤(2)数据预处理。

fMRI数据预处理：真实fMRI数据使用DPABI软件包来进行数据预处理操作。具体步骤包括：1)删除前10个时间点数据，降低实验采集过程中带来的误差。2)时间层校正，降低不同层数据因不同时间获取而造成的差异性。3)头动校正，降低被试头动带了影响。4)空间标准化，将结构项和功能像数据在被试个体水平上进行空间配准。5)数据重采样，高斯平滑。

步骤(3)选取感兴趣区域。

对于仿真数据，本发明所使用的感兴趣区域(Region of interest，ROI)共包含10个区域，10个区域的选择来源于SimTB工具包所提供的30个脑区模板。脑区的具体位置和大小如图3所示。

对于真实数据，本发明选取的感兴趣区域使用最常用的自动解剖标记(AAL)板获得大脑的90个皮层及皮层下脑区(不含26个小脑脑区)。为了探究AD被试脑效应连接网络的演化过程，直观地看到网络结构的演化，选择了116个脑区中与AD高度相关的42个脑区，脑区具体参见表1.对于每个感兴趣区域的体素时间序列，采用区域体素平均的方法，获得该区域的代表性体素。具体做法为，提取42个脑区中每个脑区(ROI)内所有体素的时间序列，并做平均，获得每个脑区的平均体素时间序列。表1 42个感兴趣区域名称

(步骤4)数据离散化。数据离散化是指将连续的fMRI数据进行分段，使其变为一段段离散化的区间。本发明采用最传统的等频离散化方法，将原始fMRI体素时间序列离散为若干个不同水平的离散化数值。等频率离散方法法即是将属性值分为具有相同宽度的区间，具体的方法为：将数据离散为d＝5部分，则每个部分被表示为低(取值为0)，较低(取值为1)，中(取值为2)，较高(取值为3)和高(取值为4)。对数据进行离散化后，数据将从原本的连续实数空间变为0-4的离散空间。

(步骤5)使用非平稳动态贝叶斯网构建与fMRI数据最匹配的动态人脑效应连接网络，具体包括以下几个步骤：

(步骤5.1)初始化参数,包括脑效应连接状态个数s，效应连接转移时间集T，每个脑区最大父母节点个数m，脑区个数n和效应连接马尔科夫阶数(延迟)l；本发明中的仿真数据参数配置如下：s＝3，T＝{60,120}，m＝5，n＝10，l＝1.真实数据的参数配置如下：s＝4，T＝{150,300,450}，m＝30，n＝42，l＝1.

(步骤5.2)构建与数据匹配的脑效应连接网络结构可以表示为：

其中G是脑效应连接网络，D是fMRI数据，P(G)是结构先验，P(D)为常数，只与fMRI数据集本身有关，对所有脑效应连接网络结构而言是一样。因此重要的是确定边缘似然函数P(D|G)。本发明采用贝叶斯评分对网络结构进行评价来计算边缘似然函数，具体表示为：

其中n为脑区个数，d为离散数量，q_i为节点脑区i的父母节点取值个数。N_ijk离散的体素时间序列x_i取值为k时，父母节点脑区取值为j的个数。

Γ(·)为伽马(Gamma)函数，α_ij和α_ijk是狄利克雷超参数先验。

将上述情况扩展到具有s个状态的脑效应连接时，可以得到：

其中G₁表示第1个状态时的动态脑效应连接网络结构，G_s表示第s个状态时的动态脑效应连接网络结构。T为效应连接转移集合，T＝{t₁,t₂,…,t_s-1}集合中的元素个数为s-1.假设从第1个状态时的动态脑效应连接网络结构G₁变为第2个状态时的G₂经历了t₁个时间点，则定义ΔG₁为网络G1到网络G2的转移网络。同样地，ΔG_s-1为第s-1个动态脑效应连接网络到第s个的动态脑效应连接网络的转移网络。引入转移网络后，公式(3)可以改写为：

P(G₁,G₂,...,G_s|T)＝P(G₁,ΔG₁,ΔG₂,...,ΔG_s-1|T) (4)

因此使用非平稳动态贝叶斯网构建动态脑效应连接网络结构也可以通过以下方式实现。首先，构建初始网络G₁，随后确定s-1个转移网络。构建过程如图4所示。

在本发明中，将脑效应连接状态个数s和效应连接时间转移集合T作为已知输入，因此公式(4)可以计算为：

其中1-e^-λs为几何先验，s_i为转移网络ΔG_i中的脑效应连接数量(有向弧数)，s_max为最大效应连接数量，

由于整个网络被分为初始效应网络和转移效应连接网络，则公式(2)，效应连接网络的BDe评分计算为：

其中h为效应连接网络的时间间隔I_h的索引值，p_i为间隔数量,q_ih为在间隔为h时，父母脑区的取值个数。从公式5中可以看出，原本父母脑区的统计数值从全部时间序列，变为间隔h的时间序列。因此，新的评分函数可以将非平稳fMRI时间划分为几部分，然后构建几部分的转移网络，最终确定整个网络结构。给定公式5的计算方式，可以通过马尔科夫链蒙特卡洛方法搜索方法(Markov chain Monte Carlo，MCMC)获得初始脑效应连接网络和转移脑效应连接网络，则动态脑效应连接被接受的概率为：

其中BDe(·)为BDe评分(公式6)，N(G)表示脑效应连接网络G邻域结构的大小，G^old表示脑效应连接网络转移(变化)前的网络，G^new表示脑效应连接网络转移(变化)后的网络。

(步骤6)结果分析。对于学得的人脑效应连接网络结构(有向无环图)，节点表示脑区(ROI)，有向弧为脑区之间的效应连接。将本发明方法与传统的动态贝叶斯网方法和贝叶斯网方法进行比较，探究新方法优势及劣势。在真实数据中，观察AD被试的动态脑效应连接网络随年龄增大而进行的变化，通过分析网络模式差异，分析AD被试脑效应连接网络的演化与年龄之间的关系。

(步骤6.1)评价指标。对比结果实验中，使用精度，准确率，召回率，F度量等四项指标，计算方式分别为：

其中TP表示算法学习到的脑效应连接与标准网络中相同效应连接的数量，TN表示算法未学习到效应连接且标准网络中也不存在该效应连接的数量，FP表示算法学习到的脑效应连接而标准网络中不存在该效应连接的数量，FN表示算法未学习到的脑效应连接而标准网络中存在该效应连接的数量。

(步骤6.2)仿真fMRI数据集上的结果与分析。

实验中，对比了5种不同的算法，它们分别为：

①PC(Peter and clerk，一种条件独立性测试贝叶斯网方法)；

②GES(Greedy Equivalence Search，基于贪婪等价类搜索的贝叶斯网方法)；

③AIAEC(Artificial immunity algorithm learning effectiveconnectivity，蚁群优化的脑效应连接网络学习方法)；

④DBN(Dynamic Bayesiannetwork，动态贝叶斯网方法)；

⑤nDBNEC(本发明提出的新方法，一种基于非平稳动态贝叶斯网的人脑效应连接网络构建方法)；

由于实验设计共包含3个状态，常规的贝叶斯网方法和DBN方法无法区分这3种不同的状态，因此将3种状态分开，分别检测4种算法在3个状态下的效果。对于本发明提出的nDBNEC方法，直接将算法在数据集上运行。实验结果如表2所示。

表2 5种算法在仿真fMRI数据集中3种状态下的实验结果对比

从表1的结果中，可以发现，传统贝叶斯网方法因为无法区分非平稳fMRI时间序列中的状态，且无法捕获fMRI数据的时间信息，因此表现不佳。动态贝叶斯网方法可以捕获时间信息，因此表现优于传统贝叶斯网方法。但是动态贝叶斯网方法的网络结构无法随着时间序列中状态的变化而变化，因此效果也受到了局限。在这种情况下，本发明提出的非平稳动态贝叶斯网方法，可以有效地捕捉fMRI数据中的非平稳时间信息，网络结构随着时间序列状态的变化而变化，因此获得了最优的结果。具体地，本发明提出的nDBNEC方法在各项指标上均显著由于其他方法，特别是对于中间状态(状态2)，其他4算法均表现不佳，而nDBNEC方法却表现出非常不错的性能。由此可见，nDBNEC方法可以有效捕捉时间信息，准确地刻画出脑效应连接网络结构随时间的变化。

(步骤6.3)真实fMRI数据集上的结果与分析

图5展示了AD被试的动态脑效应连接网络连接图，网络的结构随着年龄的增长具有4个不同的状态。具体地，t1为AD被试87岁时的脑效应连接网络，t2为该被试87.5岁时的脑效应连接网络，t3为该被试88岁时的脑效应连接网络，t4为该被试89岁时的脑效应连接网络。从动态人脑效应连接网络中，可以发现，随着年龄的增长，AD被试的脑效应连接总数有越来越少的趋势。具体来看，在被试87岁时，脑效应连接数量为105，而被试到了89岁时，脑效应连接数量减少到89。此外，一些相关脑区，颞叶(Temporal)和海马(Hippocampus)的脑效应连接，AD被试明显少于正常人。这是由于AD会引起颞叶和海马出现萎缩，记忆功能衰退。

Claims (4)

1.一种基于非平稳动态贝叶斯网的动态人脑效应连接识别方法，其特征在于：包括以下步骤，

步骤1数据获取：本方法所使用的数据，包含仿真fMRI数据和真实fMRI数据；仿真fMRI数据用开源工具包SimTB进行生成；真实fMRI数据使用ADNI数据库中的fMRI数据；

步骤2数据预处理：真实fMRI数据使用DPABI软件包来进行数据预处理；

步骤3选取感兴趣区域ROI：针对仿真fMRI数据，从DPABI软件组件模板中选择10个脑区作为感兴趣区域；对于真实fMRI数据采用AAL模板，然后提取42个与AD相关脑区的体素时间序列；

步骤4数据离散化；对连续的fMRI时间序列进行离散化处理；

步骤5使用非平稳动态贝叶斯网方法进行动态人脑效应连接网络构建；具体包括以下几个步骤：初始化参数，包括状态个数s，转移时间集t，最大父母节点个数m和马尔科夫阶数l；初始化原始时刻网络，使用马尔科夫链蒙特卡洛方法搜索初始网络中贝叶斯评分最高的贝叶斯网络；确定转移集，寻找下一时刻效应连接发生变化的节点；网络动态转换，将初始网络与转移集进行结合，确定下一时刻的贝叶斯网；重复进行网络转移过程，直到网络转移到达最后时刻；确定最终贝叶斯评分最高的动态贝叶斯网，包含一个初始人脑效应连接网络和s-1个后续人脑效应连接网络，这s个网络共同组成一个完整的、随时间变化的动态人脑效应连接网络；

步骤6结果分析：对于构建的动态人脑效应连接网络结构包含s个动态贝叶斯网络，节点表示脑区ROI，有向弧为脑区之间的效应连接；

数据离散化是指将连续的fMRI数据进行分段，使其变为一段段离散化的区间；

使用非平稳动态贝叶斯网构建与fMRI数据最匹配的动态人脑效应连接网络，具体包括以下几个步骤：

步骤5.1)初始化参数,包括脑效应连接状态个数s，效应连接转移时间集T，每个脑区最大父母节点个数m，脑区个数n和效应连接马尔科夫阶数l；仿真数据参数配置如下：s＝3，T＝{60,120}，m＝5，n＝10，l＝1.真实数据的参数配置如下：s＝4，T＝{150,300,450}，m＝30，n＝42，l＝1；

步骤5.2)构建与数据匹配的脑效应连接网络结构表示为：

其中G是脑效应连接网络，D是fMRI数据，P(G)是结构先验，P(D)为常数；采用贝叶斯评分对网络结构进行评价来计算边缘似然函数，具体表示为：

其中n为脑区个数，d为离散数量，q_i为节点脑区i的父母节点取值个数；N_ijk离散的体素时间序列x_i取值为k时，父母节点脑区取值为j的个数；

Γ(·)为伽马函数，α_ij和α_ijk是狄利克雷超参数先验；

将上述情况扩展到具有s个状态的脑效应连接时，得到：

其中G₁表示第1个状态时的动态脑效应连接网络结构，G_s表示第s个状态时的动态脑效应连接网络结构；T为效应连接转移集合，T＝{t₁,t₂,…,t_s-1}集合中的元素个数为s-1.假设从第1个状态时的动态脑效应连接网络结构G₁变为第2个状态时的G₂经历了t₁个时间点，则定义ΔG₁为第1个状态时的动态脑效应连接网络结构G1到第2个状态时的动态脑效应连接网络结构G2的转移网络；ΔG_s-1为第s-1个动态脑效应连接网络到第s个的动态脑效应连接网络的转移网络；引入转移网络后，公式(3)改写为：

P(G₁,G₂,...,G_s|T)＝P(G₁,ΔG₁,ΔG₂,...,ΔG_s-1|T) (4)

构建初始网络即第1个状态时的动态脑效应连接网络结构G₁，随后确定s-1个转移网络；

将脑效应连接状态个数s和效应连接时间转移集合T作为已知输入，公式(4)改写为：

其中1-e^-λs为几何先验，s_i为转移网络ΔG_i中的脑效应连接数量，s_max为最大效应连接数量，

公式(2)效应连接网络的BDe评分计算为：

其中h为效应连接网络的时间间隔I_h的索引值，p_i为间隔数量,q_ih为在间隔为h时，父母脑区的取值个数；基于公式5通过马尔科夫链蒙特卡洛方法搜索方法获得初始脑效应连接网络和转移脑效应连接网络，则动态脑效应连接被接受的概率为：

其中BDe(·)为BDe评分，见公式6，N(G)表示脑效应连接网络G邻域结构的大小，G^old表示脑效应连接网络转移前的网络，G^new表示脑效应连接网络转移后的网络。

2.根据权利要求1所述的一种基于非平稳动态贝叶斯网的动态人脑效应连接识别方法，其特征在于：真实fMRI数据使用DPABI软件包来进行数据预处理操作；具体步骤包括：1)删除前10个时间点数据，降低实验采集过程中带来的误差；2)时间层校正，降低不同层数据因不同时间获取而造成的差异性；3)头动校正；4)空间标准化，将结构项和功能项数据在被试个体水平上进行空间配准；5)数据重采样，高斯平滑。

3.根据权利要求1所述的一种基于非平稳动态贝叶斯网的动态人脑效应连接识别方法，其特征在于：对于仿真数据，所使用的感兴趣区域ROI共包含10个区域，10个区域的选择来源于SimTB工具包所提供的30个脑区模板。

4.根据权利要求1所述的一种基于非平稳动态贝叶斯网的动态人脑效应连接识别方法，其特征在于：对比结果实验中，使用精度，准确率，召回率，F度量四项指标，计算方式分别为：

其中TP表示算法学习到的脑效应连接与标准网络中相同脑效应连接的数量，TN表示算法未学习到脑效应连接且标准网络中也不存在该脑效应连接的数量，FP表示算法学习到的脑效应连接而标准网络中不存在该脑效应连接的数量，FN表示算法未学习到的脑效应连接而标准网络中存在该脑效应连接的数量。

Images