CN113947157B

CN113947157B - 一种基于层次聚类和结构方程模型的动态脑效应连接网络生成方法

Info

Publication number: CN113947157B
Application number: CN202111244910.0A
Authority: CN
Inventors: 冀俊忠; 韩露; 刘金铎
Original assignee: Beijing University of Technology
Current assignee: Beijing University of Technology
Priority date: 2021-10-26
Filing date: 2021-10-26
Publication date: 2024-04-30
Anticipated expiration: 2041-10-26
Also published as: CN113947157A

Abstract

本发明公开了一种基于层次聚类和结构方程模型的动态脑效应连接网络生成方法，首先使用层次聚类来评估脑状态数和转移时刻，然后使用SEM从fMRI时间序列数据中学习DEC的网络结构。从fMRI时间序列数据中获取脑状态特征信息；使用层次聚类来评估脑状态数和转移时刻；根据状态数和转移时刻对原始fMRI时间序列数据进行划分；使用SEM从fMRI时间序列数据中学习DEC的网络结构。分析在多模态数据下构建脑效应连接网络与单模态下的差异和原因，并用真实数据进行检验。本发明通过两种算法的结合，可以自动地获得动态脑效应连接网络，相较其他算法，在仿真数据集上获得了更准确的网络结构。

Description

一种基于层次聚类和结构方程模型的动态脑效应连接网络生成方法

技术领域

本发明涉及fMRI功能磁共振成像数据的脑网络生成方法，特别是一种基于层次聚类和结构方程模型的动态脑效应连接网络生成模型。

背景技术

脑效应连接定义为一个大脑区域对另一个大脑区域施加的神经影响，其损害与一些大脑疾病有关，例如阿尔茨海默病，帕金森病，和精神分裂症。其中，脑效应连接网络是由节点和有向边组成的图模型，其中节点代表脑区，有向边描述一个脑区对另一个脑区的神经活动的因果效应，与有向边相关的连接参数代表边的连接强度。脑效应连接主要包括效应连接(Effective Connectivity,EC)和动态效应连接(Dynamic EffectiveConnectivity,DEC)两种描述。

早期的脑效应连接研究侧重于在整个时间序列尺度上进行估计。例如，格兰杰因果关系(Granger Causality,GC)方法用滞后值对时间序列中的每个值进行回归，从而得到因果关系。pwLiNGAM方法根据线性非高斯非循环模型下的似然比推断效应连接。Patel方法通过贝叶斯范式比较体素的预期联合概率和边缘概率来评估不同大脑区域之间的关系。大规模格兰杰因果关系(Large-scale Granger Causality,lsGC)方法是一种数据驱动的多变量方法，其中结合了嵌入降维方法。Two-Step方法利用了BOLD信号的非高斯特征来分析大脑的效应连接。

然而，由于真实的神经元活动存在自发的时空重组现象，越来越多的方法着眼于较短的时间尺度而非整个时间序列。由于动态脑效应连接网络描述了脑区间的效应连接随时间变化的形成和消失，包含了神经活动的时间特性，所以更加符合真实的脑网络结构。DEC方法的关键问题是确定不同脑效应连接状态的转移时刻。针对这一问题，目前主要有三种解决方法：基于统计模型的方法、滑动窗口方法或参数状态估计方法。其中，统计模型的方法首先建立一个基于隐马尔可夫模型(Hidden Markov Model,HMM)和动态因果模型(Dynamic Causal Model,DCM)的概率性、生成性的血流动力学波动模型。这种统计模型可以捕捉每个时刻的脑效应连接信息。然而，这会造成维度灾难并且不利于一段时间内大脑主要活动状态的获取。滑动窗口的方法需要人为地设置滑动窗口参数，并在每个窗口下进行脑效应连接网络的构建。在滑动窗口方法中，窗口的大小和滑动步长直接决定了动态效应连接网络的规模和质量。然而，我们很难准确地选择合适的参数。参数状态估计方法应用非平稳动态贝叶斯网络(Non-stationary Dynamic Bayesian Networks,nsDBN)从功能磁共振成像数据中学习动态脑效应连接，它可以估计大脑状态的转移时刻。但是，该方法需要将连续数据离散化，这会导致数据中的信息损失。与上述两类方法相比，这种参数状态估计方法能更准确、科学地划分fMRI时间序列，从而获取更多有用的时间特征信息。然而，参数状态估计方法依然涉及到先验参数选择的问题。

针对上述问题，本发明利用层次聚类方法获得脑状态的动态信息，并使用结构方程模型(structural equation model,SEM)进行脑效应连接网络的构建，旨在解决统计模型、滑动窗口和参数状态估计方法所带来的问题，从而获得更准确的动态效应连接网络结构。

发明内容

针对现有脑效应网络构建方法的不足。本发明提出一种新的用于从fMRI时间序列数据中构建动态脑效应连接网络的生成模型。该模型能够自动地捕捉大脑状态的转移信息，并且能够表征脑区之间具有时间特征的有向依赖关系。该生成模型可以分为两部分：首先使用层次聚类来评估脑状态数和转移时刻，然后使用SEM从fMRI时间序列数据中学习DEC的网络结构。

实现本发明的主要思路是：从fMRI时间序列数据中获取脑状态特征信息；使用层次聚类来评估脑状态数和转移时刻；根据状态数和转移时刻对原始fMRI时间序列数据进行划分；使用SEM从fMRI时间序列数据中学习DEC的网络结构。分析在多模态数据下构建脑效应连接网络与单模态下的差异和原因，并用真实数据进行检验。

一种基于层次聚类和结构方程模型的动态脑效应连接网络生成模型包括以下步骤：

(步骤1)fMRI数据获取：首先为了验证本发明方法的有效性，我们使用一种fMRI数据生成方法，生成了两组仿真fMRI数据。由于已知生成的仿真数据的标准网络结构，因此可以用于验证算法有效性。

(步骤2)数据预处理：真实fMRI数据预处理使用统计参数图软件包(SPM12)来实现。

(步骤3)选取感兴趣区域(region ofinterest,ROI)：本发明采用共包含90个脑区的AAL模板，并提取这90个脑区的功能信息，即脑区的体素值(Voxels)。

(步骤4)组状态特征生成：本发明利用时变多元向量自回归(time-varyingmultivariate vector autoregressive,tvMVAR)模型，以线性差分方程的形式对脑区间的相互作用和随时间的变化进行建模，并采用双扩展卡尔曼滤波算法捕获时变模型参数。之后获得模型参数在组水平上的均值，即为脑效应连接的组状态特征。

(步骤5)使用层次聚类获取脑状态转移信息：我们使用一种无参的层次聚类算法来估计状态数和转移时刻。状态间的距离由状态特征向量的余弦相似性来评估。两个状态特征向量距离越近，表明这两个时刻的状态更相似。我们根据状态的最邻近关系对簇进行划分。计算出聚类中心的状态特征向量作为下一次迭代的聚类对象。在聚类过程停止后，我们得到了状态数和转移时刻。

(步骤6)使用结构方程模型生成动态脑效应连接网络：本发明根据状态数和转移时刻对原始fMRI时间序列数据进行划分，生成具有不同脑状态的时间序列，并利用SEM对DEC网络结构进行学习。结构方程模型的输入是每个被试每个状态下的时间序列数据，其输出是所有脑区之间的因果关系强度和模拟的时间序列。我们利用生成对抗网络(generative adversarial networks,GAN)对该模型进行训练。

结构方程模型作为生成器，判别器判断输入数据是否真实，并反馈给生成器。经过对生成器和判别的训练，我们获得了生成器产生的脑区间的效应连接。最后，当每个被试的各个状态完成训练时，该模型将得到DEC网络结构。

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

(1)本发明提出了一种新的无参状态估计方法，可以自动地估计脑状态的转移信息，获得准确的脑状态数和转移时刻。

(2)本发明应用结构方程模型进行脑效应连接网络结构学习，并利用GAN对因果关系参数进行训练，能够获得更真实、更准确的脑效应连接网络。

(3)本发明通过两种算法的结合，可以自动地获得动态脑效应连接网络，相较其他算法，在仿真数据集上获得了更准确的网络结构。

(4)本发明在真实fMRI数据集上学习的脑效应连接网络结构准确可靠，具有实际的生物学意义，可为异常效应连接的发现和脑疾病的辅助诊断提供有力的帮助。

附图说明

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

图2算法的框架。

图3时变多元向量自回归模型的框架。

图4仿真fMRI数据集的ROI示意图。

具体实施方式

下面阐述本发明的具体实施方式和详细步骤，本发明具体实现的框架如图2所示，具体包括：

(步骤1)fMRI数据获取。

为了验证本发明方法的有效性，我们使用一种常用的动态因果模型(DynamicCausal Modelling，DCM)生成了两组仿真fMRI数据。

本发明使用的真实fMRI数据来自于孤独症谱系障碍脑成像数据交换(ABIDEΙ)项目。ABIDEΙ涉及17个国际站点，共享收集的静息态功能磁共振成像。参与资料收集的患者为孤独症谱系障碍患者，对照组为健康被试。我们从纽约大学站点选取150名被试进行分类实验。患者和正常人分别占一半。

步骤(2)数据预处理。

fMRI数据的预处理操作使用统计参数图SPM软件(SPM12,http://www.fil.ion.ucl.ac.uk/spm)完成，其中包括：头动校正、时间层校正、对齐(平动和转动)、空间标准化、高斯平滑等。

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

本发明所采用的仿真fMRI数据集使用的ROI共有5个区域，如图4所示，但这5个区域并不具有实际的意义。对于真实fMRI数据，本发明选取的感兴趣区域使用最常用的自动解剖标记(AAL)模板获得大脑的90个皮层和皮层下脑区，不包含26个小脑脑区。通过对90个脑区中每个脑区(ROI)的所有体素的时间序列做平均化处理，获得每个脑区的平均BOLD信号。

(步骤4)组状态特征生成：对于时间序列从时间k的P大脑区域测量，k＝1,…,T，阶数为P的tvMVAR模型定义为：

其中[v₁(k-r)…v_P(k-r)]是白噪声向量。矩阵A_r(k)由下式给出：

对于k＝1,…,P，模型参数反映了脑区i和j在时滞参数r下的时变相互作用。

模型阶数p的值没有黄金标准。对于fMRI时间序列数据，一般认为当前时k的脑区活动仅与前一时k-1有关，因此我们将p的值设为1。我们使用双扩展卡尔曼滤波器(DEKF)来估计tvMVAR模型参数A。其等效的状态空间可以表示为：

a(k)＝a(k-1)+v_a(k) (3)

x(k)＝f(x(k-1),a(k-1))+Bv_x(k) (4)

y(k)＝Cx(k)+v_y(k) (5)

其中，我们假设待估计的模型参数a(k)由时刻k处具有额外噪声v_a(k)的非线性方程组生成。x(k)对应于由过程噪声v_x(k)驱动的真实信号，B和C分别是过程矩阵和度量矩阵。测量得到fMRI观测值是y(k)，它包含观测噪声v_y(k)。我们首先假设噪声项v_x(k)和v_y(k)是白噪声，其协方差分别是σ_x和σ_y。该估计器同时产生预测的fMRI时间序列x和模型参数a的近似最大似然估计。

我们采用两个扩展卡尔曼滤波器(EKF)，一个用于X估计，另一个用于A估计。每个EKF包含两个主要过程：时间更新和测量更新。因此，我们定义EKF₁负责X估计，产生x的预测值；EKF₂负责A估计，产生a的估计值。EKF₁与EKF₂同理，估计过程与内部结构如图3所示。EKF₁中的时间更新过程表示为：

其中是在时刻k-1下的估计向量，它乘以估计得出的以获得后验/> 是x的误差协方差矩阵。EKF₁中的测量更新过程表示为：

其中是在时刻k下的实际fMRI观测值，输出的估计值/>与/>成线性关系，系数为C。在每个时间步骤中，EKF₁计算x的卡尔曼增益/>并更新估计向量x及其误差协方差。EKF₂用于A估计，它产生A的估计值。EKF₂中的时间更新过程表示为：

其中函数vector(·)表示将矩阵展开为一维向量λ是一个超参数，通常设置为1。/>是估计向量a的误差协方差矩阵。EKF₂中的测量更新过程表示为：

其中是模型参数a的卡尔曼增益。而/>是一个动态导数，用于解释X估计滤波器的迭代过程。

对于被试i，即为第k个时刻下的组状态特征矩阵Fⁱ(k)，这将作为下一步状态划分过程的依据特征。我们使用以下公式来获得组状态特征矩阵：

对于任意时刻k,我们将矩阵F^g(k)拉伸到一维行向量f^g(k)，矩阵f^g的每一行代表每个时刻的组状态特征。

(步骤5)使用层次聚类获取脑状态转移信息：为了估计状态个数和转换时刻，我们采用了一种无阈值或其他超参数的层次聚类算法。首先，我们选择每个状态的最近邻状态，即通过计算成对余弦相似度得到的距离最小的邻居。距离矩阵$\textbf{D}$计算公式如下：

特别地，对于大规模数据，我们可以通过快速近似最近邻算法(如k-d树)获得第一邻域，而不需要计算完全成对距离矩阵。

接下来，给定每个状态的最近邻状态的整数索引，我们直接定义邻接矩阵$\textbf{M}$为：

其中，表示状态i的最近邻状态。邻接矩阵通过/>将每个状态i链接到其最近邻状态，通过/>使得矩阵强制对称，并将具有相同最近邻/>的两个状态(i,j)链接。公式(17)返回一个对称稀疏矩阵，它直接对应一个有向图，其中存在强连接的结点集是一个簇。我们将簇中心作为新的状态进行聚类，簇中心是簇中所有状态的平均值。我们可以看到，聚类方程(17)直接决定了簇的划分规则，而不依赖于任何阈值或进一步的分析。

(步骤6)使用结构方程模型生成动态脑效应连接网络：为了推断动态效应连接网络，我们采用了结构方程模型并利用由生成器和判别器组成的GAN来进行训练。在此之前，我们根据状态估计得到的状态数和转换时间对原始fMRI时间序列数据进行划分。

每个被试在每个状态下的fMRI时间序列数据作为生成器的输入。将真实的fMRI时间序列数据与生成器产生的样本配对输入GAN的判别器。判别器判断输入数据是否为真实数据，并将结果反馈给生成器。其中SEM作为生成器，用于估计脑区之间的因果关系。SEM的输入数据是划分后的fMRI时间序列数据。

其中，Θ(i)表示脑区i的父节点集，E_pi是效应连接的强度，a_pi是脑区$i$和$p$之间效应连接的因果强度参数，其方向是p到i，因此a是效应连接强度矩阵。换句话说，时间序列可以由其父代(Θ(i))时间序列和a_pi产生。另外，由于脑区和它本身之间的相互作用是没有意义的，所以我们定义：当i＝＝p时，E_pi＝0。

判别器的网络结构由一层RNN和一个全连接层(无SEM模型)组成。判别器的任务是区分生成器产生的样本和真实样本。判别器的输出表示样本是真实数据的概率。

另外，我们在效应连接生成器中添加了稀疏惩罚项，GAN的损失函数表示为：

稀疏项定义为：

其中，L是时间序列的长度，λ是控制稀疏性的超参数，||G||表示网络复杂度，计算如下：

当生成器生成的样本与实际输入数据非常相似时，我们可以通过SEM部分的因果强度参数得到效应连接。当所有状态的fMRI时间序列通过SEM训练结束后，即获得动态脑效应连接网络结构。

为了验证该模型的有效性，我们在仿真fMRI数据上与其他五种生成方法进行比较。在给定本发明估计的状态数和转换时间的情况下，我们对仿真数据集进行了划分，并对这六种EC学习方法的精度、召回率和结构汉明距离(structural hamming distance,SHD)进行了评价。

表1仿真fMRI数据集上的实验

从上表我们可以看到，我们的模型达到了最好的性能，在三种评价指标上均显著高于其他五种对比方法，其中最高精度可达到0.80。为了验证该模型在真实数据上的有效性，我们将不同的学习方法所产生的效应连接网络拉伸为特征向量，利用支持向量机(support vector machine,SVM)对患者和正常被试进行分类。我们通过分类性能的对比间接地评价效应连接方法的学习效果。我们在ABIDE I数据集上进行了三种动态效应脑连接生成方法的分类实验。

表2 ABIDE I数据集上的分类实验

方法

ACC(％)

PPV(％)

TPR(％)

TNR(％)

F1(％)

AUC(％)

滑动窗口GC

64.00±3.27

69.11±6.73

64.50±4.61

63.61±4.01

66.47±3.92

63.39±8.94

参数估计nsDBN

66.00±4.67

65.36±6.09

66.87±7.50

65.77±5.09

65.70±4.33

66.43±8.63

本发明

71.33±3.06

71.25±3.24

70.65±6.33

72.38±2.54

70.71±2.97

74.64±5.46

从表2我们可以发现，本发明在ACC、PPV、TPR、TNR、F1和AUC方面均优于对比方法。实验结果表明，本发明生成的动态效应连接网络比其他动态效应连接学习方法更能区分健康对照组和脑疾病患者。

Claims

1.一种基于层次聚类和结构方程模型的动态脑效应连接网络生成方法，其特征在于：包括以下步骤：

步骤1：fMRI数据获取：使用fMRI数据生成方法，生成两组仿真fMRI数据；

步骤2：数据预处理：真实fMRI数据预处理使用统计参数图软件包SPM12来实现；

步骤3：选取感兴趣区域ROI：采用共包含90个脑区的AAL模板，并提取这90个脑区的功能信息，即脑区的体素值Voxels；

步骤4：组状态特征生成：利用时变多元向量自回归tvMVAR模型，以线性差分方程的形式对脑区间的相互作用和随时间的变化进行建模，并采用双扩展卡尔曼滤波算法捕获时变模型参数；之后获得模型参数在组水平上的均值，即为脑效应连接的组状态特征；

步骤5：使用层次聚类获取脑状态转移信息：使用一种无参的层次聚类算法来估计状态数和转移时刻；

步骤6：使用结构方程模型生成动态脑效应连接网络：根据状态数和转移时刻对原始fMRI时间序列数据进行划分，生成具有不同脑状态的时间序列，并利用SEM对DEC网络结构进行学习；利用生成对抗网络GAN对该模型进行训练；结构方程模型作为生成器，判别器判断输入数据是否真实，并反馈给生成器；经过对生成器和判别的训练，获得生成器产生的脑区间效应连接；当每个被试的各个状态完成训练时，得到DEC网络结构。

2.根据权利要求1所述的一种基于层次聚类和结构方程模型的动态脑效应连接网络生成方法，其特征在于：步骤5中，状态间的距离由状态特征向量的余弦相似性来评估；两个状态特征向量距离越近，表明这两个时刻的状态更相似；根据状态的最邻近关系对簇进行划分；计算出聚类中心的状态特征向量作为下一次迭代的聚类对象；在聚类过程停止后，得到状态数和转移时刻。

3.根据权利要求1所述的一种基于层次聚类和结构方程模型的动态脑效应连接网络生成方法，其特征在于：步骤6中，结构方程模型的输入是每个被试者每个状态下的时间序列数据，输出是所有脑区之间的因果关系强度和模拟的时间序列。