CN105249964B - 基于脑磁图和弥散张量成像的多模态脑功能重建评估方法 - Google Patents

Abstract

本发明公开了一种基于脑磁图和弥散张量成像的多模态脑功能重建评估方法，包括如下步骤：1)在被试中选择感兴趣脑区，提取对应脑区的脑磁图信号，并对提取出的每个脑区的时间序列信号进行滑动时间窗的分窗处理，计算被试存在的功能连接矩阵；2)按照选定的感兴趣脑区，提取弥散张量成像数据，计算被试存在的结构连接矩阵；3)计算被试的功能‑结构耦合值集合；4)进行功能状态转变时间段的计算，得到最终功能连接上升速度最快的时间段；5)根据步骤3)中得到的功能‑结构耦合值集合，提取步骤4)中得到的时间段对应时间的功能‑结构耦合值。本发明方法避免了主观因素带来的误差；结合了大脑结构与功能两方面的信息，避免了单纯从功能角度判断的局限性。

Description

基于脑磁图和弥散张量成像的多模态脑功能重建评估方法

技术领域

本发明涉及一种基于脑磁图和弥散张量成像数据的多模态信号分析技术，属于图像处理和信号处理等领域。

背景技术

目前，抑郁症患者经过一段时间的治疗后，功能恢复状况的临床诊断是临床医生通过症状的缓解，结合汉密尔顿抑郁量表(HAMD)、大体评定量表(GAS)等量表，通过医师评估及患者自评得到的。这种方式，具有很强的主观性，对临床医生的问询经验有极高的要求。另外，症状的外显往往会迟于脑环路的改变。这些因素将导致难以客观及时地反映患者的真实恢复程度。

抑郁症患者经过急性期治疗，功能活动水平可以得到一定程度的恢复，但是可能会出现反复发作情况，即功能并没有得到真正的恢复。既往对抑郁症患者功能恢复的研究，多是单独从功能角度来研究，然而大脑的结构与功能是紧密联系的，功能是建立在大脑结构的基础上，结构与功能之间关系的异常可能提示着抑郁症功能异常的物质基础。因此，从脑结构和功能影像学中提取多模态量化指标，可以为治疗后大脑功能活动水平的恢复程度，提供更为有效和及时的评估，从而保证足程足量的有效治疗，同时避免过度治疗。

发明内容

发明目的：针对现有技术中存在的问题，本发明提出一种基于脑磁图和弥散张量成像的多模态脑功能重建评估方法，从功能、结构两方面融合的角度及时有效评估抑郁症患者功能重建情况，避免主观因素干扰。

技术方案：为了解决上述问题，本发明提出一种基于脑磁图和弥散张量成像的多模态脑功能重建评估方法，首先依靠脑磁图数据获得动态脑功能连接信息，依靠弥散张量成像数据获得脑结构连接信息；然后对感兴趣脑区间的脑功能连接时间序列进行分析，获取脑区间功能连接发生快速上升的阶段，在该阶段意味着为了进行有效信息交互，功能连接需要摆脱结构的约束来更加柔性地处理当前任务；最后计算功能连接发生快速上升的阶段的功能结构耦合，以此标记脑功能为了有效处理任务的柔性能力，同时也反映患者通过治疗获得的功能重建程度，具体包括如下步骤：

1)在被试中选择感兴趣脑区，提取对应脑区的脑磁图信号，并对提取出的相应信号的时间序列进行滑动时间窗的分窗处理，计算对应时间窗内每两个脑区间的功能连接值，进而得到被试存在的功能连接矩阵，具体计算方法为：

在每个时间窗内，对每个脑区内的信号进行平均化，利用皮尔逊相关计算每两个脑区之间的时间相关性，作为该时间窗内这两个脑区间的功能连接值，进而得到这两个脑区间的功能连接矩阵。每个被试存在一个M×M×N的功能连接矩阵，其中M是感兴趣脑区的数目，N是滑动时间窗的个数，每个时间窗的窗长为：T_j＝[t₁+(j-1)×t₄,t₁+(j-1)×t₄+t₃],j＝1,2,…,t₁,t₂,t₃,t₄分别为时间序列的起始时间、结束时间、时间窗长度和滑动步长。

2)按照步骤1)中选定的感兴趣脑区，提取弥散张量成像数据，计算每两个脑区间的结构连接值，进而得到被试存在的结构连接矩阵，具体计算方法为：

对于任意两个脑区Roi(u)和Roi(v)，可视为节点u和v，其皮层面积分别为S(u)和S(v)，在结构网络中连接节点u和v的边e(u,v)对应于连接Roi(u)和Roi(v)的神经纤维，边的权重为：

其中，f为节点u和v之间的一条纤维，E_f为连接节点u和v的所有纤维，l(f)为这些纤维的长度，w(e)反映的是两区域的连接密度，以此作为结构连接值，因此每个被试存在一个M×M的结构连接矩阵。为了与步骤1)中建立的功能连接矩阵相匹配，将结构连接矩阵重复N次，得到结构连接矩阵M×M×N，N是滑动时间窗的个数，每个窗口下的结构连接矩阵相同。

3)分别计算被试在每个时间窗内的功能-结构耦合值，进而得到被试的功能-结构耦合值集合，具体计算方法为：

将被试对应的功能连接矩阵和结构连接矩阵转换为两列向量，计算两列向量的皮尔逊相关值，即为该时间窗内的功能-结构耦合值，最终得到功能-结构耦合值集合{fs_i},i＝1,2,…,N，N是滑动时间窗的个数。

4)进行功能状态转变时间段的计算，得到最终功能连接上升速度最快的时间段，具体求解步骤为：

4.1)提取步骤1)中任意两个脑区间的功能连接序列将功能连接序列分为[1,γ]和[γ+1,N]两部分，分别计算这两部分及整体的贝叶斯信息量BIC，BIC＝-ln(L)+ln(n)×k，其中L是最大似然估计，n是计算贝叶斯信息量的功能连接序列的数据点个数，k是对应的功能连接序列的参数个数，M是感兴趣脑区的个数，N是滑动时间窗的个数；重复上述步骤从γ＝2到γ＝N-1，求解对应的γ’，即为可能状态转变时间点，此时功能连接序列被分为[1,γ’]和[γ’+1,N]两部分，其中BIC_1γ,BIC_2γ,BIC分别为[1,γ]，[γ+1,N]和[1,N]的贝叶斯信息量；分别对[1,γ’]和[γ’+1,N]重复上述步骤，直到无法继续分割为止，得到可能状态转变时间点集合{γ'_k1},k1＝1,2,…,s₁，s₁是可能状态转变时间点个数；

4.2)将功能连接序列[1,γ’_k1]分为[1,γ’_k1-1]和[γ’_k1-1+1,γ’_k1]两部分，若这两部分的贝叶斯信息量之和小于功能连接序列[1,γ’_k1]的贝叶斯信息量，则保留可能状态转变时间点γ’_k1-1，k1＝2,…,s₁，否则舍去；

4.3)重复步骤4.2)，直到不再有可能状态转变时间点被去除，得到新的可能状态转变时间点集合{γ'_k22},k2＝1,2,…,s₂，s₂是新的可能状态转变时间点个数；

4.4)将[1,N]时间点上的功能连接序列平均分为R个分块，每个分块即为一个功能连接变化分块，R是L的正约数。利用bootstrap检验得到具有显著性的状态转变时间点集合{γ'_k3},k3＝1,2,…,s₃，s₃是确定的显著性状态转变时间点个数；

4.5)分别计算每一个包含显著性状态转变时间点γ'_k3的功能连接变化分块的一阶导数之和，提取一阶导数之和最大的分块T_max，T_max表示功能连接上升速度最快的时间段；

4.6)对于所有脑区对之间的功能连接序列，重复步骤4.1)到4.5)，得到集合M为感兴趣脑区个数，{T_jmax}表示功能连接上升速度最快的时间段集合，取集合{T_jmax}的众数T_mode做为最终功能连接上升速度最快的时间段。

5)根据步骤3)中得到的功能-结构耦合值集合，提取步骤4)中得到的时间段对应时间的功能-结构耦合值，该值便是评价该被试的功能重建的多模态量化指标。

有益效果：本发明采用如上技术方案，具有以下优点：

1)纯数据驱动，完全由抑郁症患者脑区信号研究来判断患者的功能恢复程度，此过程没有人工判断参与，避免了主观因素带来的误差。

2)结合了大脑结构与功能两方面的信息，从更深层次判断判断患者的功能恢复程度，避免了单纯从功能角度判断的局限性。

3)为经过急性期治疗，功能活动水平得到恢复的患者提供了疗效判断量化手段，为临床优化治疗提供辅助信息。

附图说明

图1为本发明的方法流程示意图；

图2为本发明实施例在功能连接上升速度最快时间段的功能连接和功能-结构耦合值的组间比较图，其中图2-a为治疗前患者组、治疗后患者组及健康对照组在功能连接上升速度最快时间段的功能连接值比较图，图2-b为治疗前患者组、治疗后患者组及健康对照组在功能连接上升速度最快时间段的功能-结构耦合值比较图，图中线段的上下界为功能连接和功能-结构耦合值的最大值和最小值，‘*’、‘o’、‘□’分别为三组的平均值，MDD1＝治疗前患者组，MDD2＝治疗后患者组，HC＝健康对照组；

图3为本发明实施例的治疗后患者的功能重建评估图，图中MDD1＝治疗前患者组，MDD2＝治疗后患者组，HC＝健康对照组，‘*’、‘○’、‘□’分别为三组的平均值，‘☆’表示进行功能重建评估的治疗后患者。

具体实施方式

下面结合具体实施例，进一步阐明本发明，应理解这些实施例仅用于说明本发明而不用于限制本发明的范围，在阅读了本发明之后，本领域技术人员对本发明的各种等效变换均落于本申请所附权利要求所限定的范围。

本实施例选择26例未经治疗的抑郁症患者、16例治疗后的抑郁症患者，26例性别、年龄、受教育程度相匹配的健康对照者，三组受试进行脑磁图和弥散张量成像信号采集，并进行预处理及源重建。选取额顶网络为研究对象，包括左侧顶内沟、右侧顶内沟、左侧前额叶、右侧前额叶、左侧楔前叶、右侧楔前叶、中扣带回、左侧顶下小叶、右侧顶下小叶、左背外侧前额叶和右背外侧前额叶。提取0-900ms的时间序列，以100毫秒为窗长，10毫秒为步长进行动态分窗，分别计算每个时间窗内的功能连接和结构连接，然后计算得到功能-结构耦合值。计算功能状态转变时间段，并提取对应时间段的功能连接值和功能-结构耦合值。

图2为本实施例在功能连接上升速度最快时间段的功能连接和功能-结构耦合值的组间比较图，由图2a可知，治疗后患者组的功能连接水平接近正常对照，显著高于未治疗患者组，表明治疗后功能活动水平得到好转；但图2b治疗后患者组的功能-结构耦合值并未有显著改善，表明功能与结构的关联性仍存在异常，需进一步治疗。

图3为本实施例的治疗后患者的功能重建评估图，显示了我们所提出的客观评价量化指标在临床中的使用示意，根据图示，我们通过新测试个体的指标的位置来估计被试的疾病恢复情况。由图3可以看出，该患者的功能-结构耦合值较健康对照偏高，与治疗前患者组相比无明显差异，需要进一步治疗。

Claims (5)

1.一种基于脑磁图和弥散张量成像的多模态脑功能重建评估方法，其特征在于，包括如下步骤：

1)在被试中选择感兴趣脑区，提取对应脑区的脑磁图信号，并对提取出的相应信号的时间序列进行滑动时间窗的分窗处理，计算对应时间窗内每两个脑区间的功能连接值，进而得到被试存在的功能连接矩阵；

2)按照步骤1)中选定的感兴趣脑区，提取弥散张量成像数据，计算每两个脑区间的结构连接值，进而得到被试存在的结构连接矩阵；

3)分别计算被试在每个时间窗内的功能-结构耦合值，进而得到被试的功能-结构耦合值集合；

4)进行功能状态转变时间段的计算，得到最终功能连接上升速度最快的时间段；

5)根据步骤3)中得到的功能-结构耦合值集合，提取步骤4)中得到的时间段对应时间的功能-结构耦合值，该值为评价该被试的功能重建的多模态量化指标。

2.如权利要求1所述的基于脑磁图和弥散张量成像的多模态脑功能重建评估方法，其特征在于，所述对应时间窗内每两个脑区间的功能连接值的计算方法为：在每个时间窗内，对每个脑区内的信号进行平均化，利用皮尔逊相关计算每两个脑区之间的时间相关性作为该时间窗内这两个脑区间的功能连接值。

3.如权利要求1所述的基于脑磁图和弥散张量成像的多模态脑功能重建评估方法，其特征在于，所述每两个脑区间的结构连接值的计算方法为：

对于任意两个脑区Roi(u)和Roi(v)，可视为节点u和v，其皮层面积分别为S(u)和S(v)，在结构网络中连接节点u和v的边e(u,v)对应于连接Roi(u)和Roi(v)的神经纤维，这两个脑区间的结构连接值为边的权重，即：

<mrow> <mi>w</mi> <mrow> <mo>(</mo> <mi>e</mi> <mo>)</mo> </mrow> <mo>=</mo> <mfrac> <mn>2</mn> <mrow> <mi>S</mi> <mrow> <mo>(</mo> <mi>u</mi> <mo>)</mo> </mrow> <mo>+</mo> <mi>S</mi> <mrow> <mo>(</mo> <mi>v</mi> <mo>)</mo> </mrow> </mrow> </mfrac> <msub> <mi>&amp;Sigma;</mi> <mrow> <mi>f</mi> <mo>&amp;Element;</mo> <msub> <mi>E</mi> <mi>f</mi> </msub> </mrow> </msub> <mfrac> <mn>1</mn> <mrow> <mi>l</mi> <mrow> <mo>(</mo> <mi>f</mi> <mo>)</mo> </mrow> </mrow> </mfrac> </mrow>

其中，f为节点u和v之间的一条纤维，E_f为连接节点u和v的所有纤维，l(f)为这些纤维的长度。

4.如权利要求1所述的基于脑磁图和弥散张量成像的多模态脑功能重建评估方法，其特征在于，所述被试在每个时间窗内功能-结构耦合值的计算方法为：将每个被试对应的功能连接矩阵和结构连接矩阵转换为两列向量，计算两列向量的皮尔逊相关值，即为该时间窗内的功能-结构耦合值。

5.如权利要求1所述的基于脑磁图和弥散张量成像的多模态脑功能重建评估方法，其特征在于，所述最终功能连接上升速度最快的时间段的具体求解步骤为：

4.1)提取步骤1)中任意两个脑区间的功能连接序列将功能连接序列分为[1,γ]和[γ+1,N]两部分，分别计算这两部分及整体的贝叶斯信息量BIC，BIC＝-ln(L)+ln(n)×k，其中L是最大似然估计，n是计算贝叶斯信息量的功能连接序列的数据点个数，k是对应的功能连接序列的参数个数，M是感兴趣脑区的个数，N是滑动时间窗的个数；重复上述步骤从γ＝2到γ＝N-1，求解对应的γ’，即为可能状态转变时间点，此时功能连接序列被分为[1,γ’]和[γ’+1,N]两部分，其中BIC_1γ,BIC_2γ,BIC分别为[1,γ]，[γ+1,N]和[1,N]的贝叶斯信息量；分别对[1,γ’]和[γ’+1,N]重复上述步骤，直到无法继续分割为止，得到可能状态转变时间点集合{γ'_k1},k1＝1,2,…,s₁，s₁是可能状态转变时间点个数；

4.2)将功能连接序列[1,γ’_k1]分为[1,γ’_k1-1]和[γ’_k1-1+1,γ’_k1]两部分，若这两部分的贝叶斯信息量之和小于功能连接序列[1,γ’_k1]的贝叶斯信息量，则保留可能状态转变时间点γ’_k1-1，k1＝2,…,s₁，否则舍去；

4.3)重复步骤4.2)，直到不再有可能状态转变时间点被去除，得到新的可能状态转变时间点集合{γ'_k2},k2＝1,2,…,s₂，s₂是新的可能状态转变时间点个数；

4.4)将[1,N]时间点上的功能连接序列平均分为R个分块，每个分块即为一个功能连接变化分块，R是L的正约数，利用bootstrap检验得到具有显著性的状态转变时间点集合{γ'_k3},k3＝1,2,…,s₃，s₃是确定的显著性状态转变时间点个数；

4.5)分别计算每一个包含显著性状态转变时间点γ'_k3的功能连接变化分块的一阶导数之和，提取一阶导数之和最大的分块T_max，T_max表示功能连接上升速度最快的时间段；

4.6)对于所有脑区对之间的功能连接序列，重复步骤4.1)到4.5)，得到集合M为感兴趣脑区个数，{T_jmax}表示功能连接上升速度最快的时间段集合，取集合{T_jmax}的众数T_mode做为最终功能连接上升速度最快的时间段。