CN102509282B - 一种融合结构连接的各脑区间的效能连接分析方法 - Google Patents

Abstract

一种融合结构连接的各脑区间的效能连接分析方法，步骤是提取并分析感兴趣区的大脑结构网络，并将得到的结构信息转换到效能连接参数的先验概率分布空间；然后建立基于变分贝叶斯框架的效能连接模型；最后通过集成学习和EM算法，求取各脑区间的效能连接。本发明相对于其他方法具有以下优点：1、通过转换模型将结构连接映射到效能连接参数的先验概率空间，并在后续的集成学习中优化模型参数，使得结构连接和效能连接的关系得到真实的反映；2、结合了结构连接信息，使得脑活动分析的结果更加可靠，而且便于对个体进行实际情况的探讨。

Description

一种融合结构连接的各脑区间的效能连接分析方法

技术领域

本发明涉及一种融合结构连接的效能连接计算方法，属于医学影像方法学和医学影像信号处理等领域，适用于大脑活动的分析和研究。

背景技术

既往研究发现包括抑郁症、精神***症、癫痫等精神疾病都与大脑的结构的受损或者功能性障碍有关。大脑的结构信息(结构连接)与功能信息(功能连接、效能连接)反映了整个大脑的状况，并且对各种精神疾病的病理机制的探索和研究起着重要的作用。而且结构连接和功(效)能连接之间存在密切的联系，功(效)能连接在一定的程度上反映了结构连接，但不完全取决于结构连接。以往的研究很少将它们结合起来进行研究，若从多模态的角度融合结构连接和功能连接，将会给各种精神疾病的研究带来极大的帮助。

关于结合结构连接和效能连接进行大脑活动的研究，已有一些初步的方法和探索但存在一些问题。方法1：利用DTI(弥散张量磁共振成像diffusion tensorimaging)和MEG(脑磁图magnetoencephalography)各自的优点，从不同的角度单独进行研究分析，然后综合分析。但是没有探讨这两者的关系，更没涉及到它们之间的融合。方法2：运用对结构连接的分析对功能连接的结果进行约束和优化。该方法虽然运用并分析了两者之间的关系，但是仍然没有能够真正将两者融合起来。方法3：利用结构连接的实验数据来大规模模拟大脑的功能，或是研究什么类型的网络结构会产生特定功能的神经活动。该方法能够将两者融合起来，但是仍只是重在探讨两者之间的联系。

发明内容

为了解决现有方法中的不足之处，本发明提出了一种新的融合结构连接的各脑区间的效能连接分析方法，具体技术方案如下：

融合结构连接的各脑区间的效能连接分析方法，基本思想是：提取并分析感兴趣区的大脑结构网络，并将得到的结构信息转换为效能连接参数的先验概率分布空间；然后建立基于变分贝叶斯框架的效能连接模型；最后通过集成学习(Ensemble Learning)和EM算法(Expectation-Maximization Algorithm)，求取各脑区间的效能连接。

本方法的具体步骤包括：

1、一种融合结构连接和效能连接进行脑活动分析的方法，其特征在于步骤包括：

1)首先利用弥散张量磁共振成像DTI数据进行全脑的神经纤维追踪，并建立整个大脑的结构连接网络；另外，对采集得到的脑磁图MEG信号进行3D源重建；MEG的3D源重建运用了MSP(Multiple Sparse Priors)的方法，该步骤主要是要SPM8(http://www.fil.ion.ucl.ac.uk/spm/)软件进行处理。该方法是一种运用分层贝叶斯及经验贝叶斯进行分布式偶极子源重建的方法，其的优点是能够在简单的先验知识下进行多个稀疏皮层源自动选取。

2)根据需要分析的感兴趣区域提取出所述连接网络中的主干网络，并将主干网络转化为图；

设主干网络中的任意感兴趣脑区为i，主干网络中的每个感兴趣区ROI(i)均视为一个节点i，节点i代表的感兴趣区域的皮层面积为S(i)；

主干网络中连接两个脑区i和j的神经纤维ROI(i)和ROI(j)对应于连接节点i和j的边E(i，j)，边的长度和权重分别为

和

其中，E_f为连接节点i和j的所有纤维，l_f为这些纤维的长度，N_f为神经纤维的数目；所以l为连接两区域的所有神经纤维的平均长度，w反映的是两区域的连接密度；

3)对所得的各个感兴趣区的结构信息(这里的结构信息指的是边的长度和权重，对应于连接感兴趣区的神经纤维的长度和密度)进行归一化，并转换到效能连接的先验概率分布空间，转换模型为： ${\mathrm{\Σ}}_{\mathrm{ij}}=\frac{{\mathrm{\Σ}}_0}{1+{\mathrm{\Σ}}_0{e}^{a-{\mathrm{bs}}_{\mathrm{ij}}}}={\mathrm{\α}}_{\mathrm{ij}}^{-1}={\mathrm{\α}}_0+e^{a-{\mathrm{bs}}_{\mathrm{ij}}};$

任意脑区i和j之间的效能连接服从高斯分布N(0，∑_ij)，s_ij为归一化后的结构连接信息，∑₀、a、b为模型的可调节参数；

4)基于变分贝叶斯框架的效能连接模型为Y＝XW+E，模型Y中，W对应表示脑区间的效能连接参数矩阵，该矩阵W为自回归系数参数矩阵；

E是均值为零、精度矩阵为Λ的高斯噪声，且Λ～Γ(b，c)；X，Y为经过3D源重建后的感兴趣区域信号，对于给定的数据集D＝{X，Y}有： $p\left(D\right|W,\mathrm{\Λ})={\left(2\mathrm{\π}\right)}^{-\mathrm{dN}/2}{\left|\mathrm{\Λ}\right|}^2{e}^{-\frac{1}{2}\mathrm{Tr}\left(\mathrm{\Λ}{E}_D\left(W\right)\right)};$ 为了便于模型Y的分析，将W拉长为向量w，w的分布如下：

$p\left(w\right|\left\{{\mathrm{\α}}_k\right\})=\underset{k=1}{\overset{n}{\mathrm{\Π}}}{\left(\frac{{\mathrm{\α}}_k}{2\mathrm{\π}}\right)}^{1/2}{e}^{-{\mathrm{\α}}_k{E}_k\left(w\right)};$ 其中，d为脑区信号的个数，N为脑区信号序列的长度，w为矩阵W拉伸的向量，n为效能连接参数的个数， $E_k\left(w\right)=\frac{1}{2}{w}^T{I}_{k}w;$

所述参数w、Λ和α服从高斯分布N(0，∑_ij)；

5)通过集成学习方法和最大期望EM算法，求取各脑区间的效能连接Y。

有益效果

本发明相对于其他方法具有以下优点：1、通过转换模型将结构连接映射到效能连接参数的先验概率空间，并在后续的集成学习中优化模型参数，使得结构连接和效能连接的关系得到真实的反映；2、结合了结构连接信息，使得脑活动分析的结果更加可靠，而且便于对个体进行实际情况的探讨。

附图说明

图1：本方法的基本流程示意图；

图2：结构连接信息处理流程示意图；

图3：结构信息转换模型图示意图；

图4：效能连接结果示意图。

具体实施方式

现结合附图对本发明作进一步的描述：

本发明的整个流程可以参考附图1，具体的实施步骤如下：

1、大脑结构连接网络的建立以及脑磁信号的预处理和源重建

从核磁共振弥散图像DTI到整个大脑的高精度结构连接网络的生成需要以下几个步骤：(1)弥散加权图像的处理，如涡流校正、头动校正，弥散张量模型的拟合、弥散张量和各向异性值的计算，以及标准空间的转换等；(2)灰白质的分割；(3)白质神经纤维束追踪成像；(4)大脑皮层结构的分割和感兴趣区域的选取；(5)根据需要分析的感兴趣区域提取出的主干网络，并将其转化为图。网络中的每个感兴趣区ROI(i)均可视为一个节点i，其代表的区域的皮层面积为S(i)。网络中连接ROI(i)和ROI(j)的神经纤维对应于连接节点i和j的边E(i，j)，边的长度和权重分别为 $l=\frac{1}{N_f}{\mathrm{\Σ}}_{f\∈E_f}{l}_f,$ $w=\frac{2}{S\left(i\right)+S\left(j\right)}{\mathrm{\Σ}}_{f\∈E_f}\frac{1}{l_f}.$ 其中，E_f为连接节点i和j的所有纤维，l_f为这些纤维的长度，N_f为神经纤维的数目，故l指的是连接两区域的所有神经纤维的平均长度，w反映的是两区域的连接密度。该处理流程可参考附图2。

脑磁信号的预处理(包括转化、分割、滤波、去除伪迹、平均化)后，进行3D源重建，该处理流程使用了MSP(Multiple Sparse Priors)的方法，该方法步骤主要是要SPM8(http://www.fil.ion.ucl.ac.uk/spm/)软件进行处理。该方法是一种运用分层贝叶斯及经验贝叶斯进行分布式偶极子源重建的方法，其的优点是能够在简单的先验知识下进行多个稀疏皮层源自动选取。

2、结构连接信息的转化

将上述所得的各个感兴趣区的结构信息进行归一化，并将其转换到效能连接的先验概率分布空间，转换模型为：

${\mathrm{\Σ}}_{\mathrm{ij}}=\frac{{\mathrm{\Σ}}_0}{1+{\mathrm{\Σ}}_0{e}^{a-{\mathrm{bs}}_{\mathrm{ij}}}}={\mathrm{\α}}_{\mathrm{ij}}^{-1}={\mathrm{\α}}_0+e^{a-{\mathrm{bs}}_{\mathrm{ij}}}$

任意脑区i和j之间的效能连接服从高斯分布N(0，∑_ij)，s_ij为归一化后的结构连接信息，∑₀、a、b为模型的可调节参数，在后续的集成学习中逐步地对它们进行优化，使得结构连接和效能连接的关系得到真实的反映。该转换模型的函数图如附图3所示。

3、基于变分贝叶斯框架的效能连接模型的分析

基于变分贝叶斯框架的效能连接模型为Y＝XW+E。模型中W为自回归系数矩阵，对应于脑区间的效能连接参数矩阵，这些参数服从高斯分布N(0，∑_ij)；E是均值为零，精度矩阵为Λ的高斯噪声，且Λ～Γ(b，c)。X，Y为经过源重建后的感兴趣区域信号，对于给定的数据集D＝{X，Y}有：

$p\left(D\right|W,\mathrm{\Λ})={\left(2\mathrm{\π}\right)}^{-\mathrm{dN}/2}{\left|\mathrm{\Λ}\right|}^2{e}^{-\frac{1}{2}\mathrm{Tr}\left(\mathrm{\Λ}{E}_D\left(W\right)\right)}$

为了便于模型的分析，将系数参数矩阵W拉长为以向量w，其分布如下：

$p\left(w\right|\left\{{\mathrm{\α}}_k\right\})=\underset{k=1}{\overset{n}{\mathrm{\Π}}}{\left(\frac{{\mathrm{\α}}_k}{2\mathrm{\π}}\right)}^{1/2}{e}^{-{\mathrm{\α}}_k{E}_k\left(w\right)}$

其中d脑区信号的个数，N为信号序列的长度，w为效能连接参数矩阵W拉伸的向量，n为效能连接参数的个数，

4、模型的求解，本例中给出一种改进的EM的算法：

对于以上给定的数据集D和参数θ＝{w，α，Λ}，模型m的对数论据为

log p(D|m)＝F(θ)+KL(q(θ|D)||p(θ|D，m))

F(θ)为模型的负自由能量，当q(θ|D)＝p(θ|D，m)，即模型参数的近似后验概率分布等同于真实后验概率分布时模型对数论据取得下限F(θ)，此时的参数也正是模型要求参数。然而直接通过等式求解参数是十分困难的，将F(θ)进一步分解得到F(θ)＝∫q(θ|D)log p(D |θ，m)dθ-KL(q(θ|D)||p(θ|m))。给定参数θ的初始值并固定其中的α和Λ，通过简化可得当q(w|D)＝e^I(w)时，F(θ)取得最大值，其中I(w)＝∫∫q(Λ|D)q(α|D)log(p(D|w，Λ)p(w|α))dαdΛ。更新参数w并固定w和Λ，同理可求解出α，然后更新参数α并固定α和w求解出Λ。通过这种改进的EM(Expectation-Maximization)的算法，不断进行迭代直至收敛，求出模型的参数和各脑区之间的效能连接(结果见附图4)。

Claims (3)

1.一种融合结构连接的各脑区间的效能连接分析方法，其特征在于步骤包括：

1）首先利用弥散张量磁共振成像DTI数据进行全脑的神经纤维追踪，并建立整个大脑的结构连接网络；另外，对采集得到的脑磁图MEG信号进行3D源重建；

2）根据需要分析的感兴趣区域提取出所述连接网络中的主干网络，并将主干网络转化为图；

设主干网络中的任意感兴趣脑区为i，主干网络中的每个感兴趣区ROI(i)均视为一个节点i，节点i代表的感兴趣区域的皮层面积为S(i)；

主干网络中连接两个脑区i和j的神经纤维ROI(i)和ROI(j)对应于连接节点i和j的边E(i,j)，边的长度和权重分别为

和

其中，E_f为连接节点i和j的所有纤维，l_f为这些纤维的长度，N_f为神经纤维的数目；所以l为连接两区域的所有神经纤维的平均长度，u反映的是两区域的连接密度；f为连接脑节点i和j的一根神经纤维；

3）对所得的各个感兴趣区的结构信息进行归一化，并转换到效能连接的先验概率分布空间，转换模型为：

所述结构信息指的是所述边的长度和权重，对应于连接感兴趣区的神经纤维的长度和密度；

任意脑区i和j之间的效能连接服从高斯分布N(0,Σ_ij)，s_ij为归一化后的结构连接信息，Σ₀、a、b为模型的可调节参数；其中，

α₀是常数参数，参数a,b控制结构连接对先验概率分布方差的约束程度，a为拐点参数，b为斜率参数；

4）基于变分贝叶斯框架的效能连接模型为Y=XW+E，模型Y中，W对应表示脑区间的效能连接参数矩阵，该矩阵W为自回归系数参数矩阵；X为时间序列矩阵；

E是均值为零、精度矩阵为Λ的高斯噪声，精度矩阵Λ对应元素Λ_ij服从如下高斯分布：Λ_ij～N(0,Σ_ij)；X，Y为经过3D源重建后的感兴趣区域信号，对于给定的数据集D={X,Y}有： $p\left(D\right|W,\mathrm{\Λ})={\left(2\mathrm{\π}\right)}^{-\mathrm{dN}/2}{\left|\mathrm{\Λ}\right|}^2{e}^{-\frac{1}{2}\mathrm{Tr}\left({\mathrm{\ΛE}}_D\left(W\right)\right)},$ E_D(Ｗ)＝(Y-XW)^T(Y-XW)为未正规化的误差协方差矩阵，|Λ|²为矩阵Λ行列式的平方；为了便于模型Y的分析，将W拉长为向量w，w的分布如下：

其中，d为脑区信号的个数，N为脑区信号序列的长度，w为矩阵W拉伸的向量，n为效能连接参数的个数，

α_k为第k个效能连接参数的参数精度；I_k为第k个单位矩阵；

所述参数w_ij,λ_ij,α_ij服从高斯分布N(0,Σ_ij)；

5）通过集成学习方法和最大期望EM算法，求取各脑区间的效能连接Y。

2.根据权利要求1所述的方法，其特征是所述步骤1）中，建立整个大脑的结构连接网络的步骤包括：

101）DTI图像进行处理，包括：

涡流校正和头动校正、弥散张量模型的拟合、弥散张量和各向异性值的计算、以及标准空间的转换；

102）灰白质的分割；

103）白质神经纤维束追踪成像；

104）大脑皮层结构的分割和感兴趣区域的选取。

3.根据权利要求1所述的方法，其特征是所述步骤5）中，对于所述数据集D和参数θ={w,α,Λ}，模型Y的对数论据为：logp(D|m)=F(θ)+KL(q(θ|D)||p(θ|D,m))，其中F(θ)为模型Y的负自由能量，当q(θ|D)=p(θ|D,m)，即模型参数的近似后验概率分布等同于真实后验概率分布时，模型Y对数论据取得下限F(θ)时，参数θ正是模型Y要求参数；其中，m为产生数据集D的模型；KL(q(θ|D)||p(θ|D,m))表示q(θ|D)与p(θ|D,m)的KL散度，KL散度是Kullback-Leibler差异（Kullback-Leibler Divergence）的简称，也叫做相对熵（Relative Entropy），其物理意义是：在相同事件空间里，概率分布q(θ|D)的事件空间，若用概率分布p(θ|D,m)编码时，平均每个基本事件符号编码长度增加了多少比特；是距离的表示符号；

将F(θ)进一步分解得到：

F(θ)=∫q(θ|D)logp(D|θ,m)dθ-KL(q(θ|D)||p(θ|m))；

给定参数θ的初始值并固定其中的α和Λ，通过简化得到：当q(w|D)=e^I(w)时，F(θ)取得最大值，其中I(w)=∫∫q(Λ|D)q(α|D)log(p(D|w,Λ)p(w|α))dαdΛ；更新参数w并固定w和Λ，同理求解出α，然后更新参数α并固定α和w，求解出Λ；

通过该算法，不断进行迭代直至收敛，求出模型的参数和各脑区之间的效能连接。

Images