一种基于功能磁共振心理生理交互的脑网络构建方法与流程

本发明涉及脑网络构建领域，尤其涉及一种基于功能磁共振心理生理交互的脑网络构建方法。

背景技术：

探索人脑与行为学间的对应关系最为流行的技术手段就是基于血氧水平依赖(bold,blood-oxygen-leveldependent)的功能磁共振成像。而对行为学任务的功能磁共振成像进行分析的方法有很多，其中大部分是计算局部脑信号在不同行为学状态下的改变从而推导出与行为相关的脑活动区。如因子设计统计模型法，计算任务态因子与静息态的bold信号差异了解与任务相关的脑活动。但通常行为学任务不是由一个脑区单独完成的，而是有很多个脑区共同协作，局部分析存在对行为学水平认识不全等问题。因此从全局角度出发分析脑区间的功能协作模式能让我们更全面地了解大脑的行为活动机制。

因此，针对局部分析存在的问题，提出了全脑网络分析方法。全脑网络分析方法在脑网络分析中被广泛使用，尤其是无任务刺激下(静息态)的大脑固有脑网络研究。但是现有的技术方案存在以下不足：对于任务刺激下的功能磁共振信号，直接计算各个脑区间的bold时间相关系数值并不能很好的反映出这些所构建的网络拓扑组织是否与任务状态相关，因为所提取到的各个脑区间的bold信号涵盖了脑区所有的激活信号，并不能确定这些bold信号与具体某个任务参数有关。所构建的脑网络与任务参数的对应关系不清楚，不能明确指出哪些网络连接与行为表现相关。

因此，如何构建一个与任务参数对应关系清楚的脑网络成为脑网络构建领域亟需解决的问题。

技术实现要素：

针对现有技术之不足，本发明提出一种基于功能磁共振心理生理交互的脑网络构建方法，所述方法包括：

步骤1：采集被试大脑在做任务状态下的磁共振数据，所述磁共振数据包括结构磁共振数据和功能磁共振数据，其中功能磁共振数据是一个四维的时间序列图；

步骤2：对所述功能磁共振数据进行预处理，步骤如下：

步骤21：去除所述功能磁共振数据前四个时间点的数据；

步骤22：时间层校正，采用辛格插值法对所述功能磁共振数据进行时间偏移量校正；

步骤23：头动校正，将时间层校正后的中间层图作为基准图像，然后通过旋转或平移等刚性变换将时间序列上的其它层图像与第一层图像的位置匹配，然后用内插值算法对磁共振数据进行重新采样，同时至少生成平均值文件和头动参数；

步骤24：减少运动伪影，若每个被试者的功能磁共振数据头动的任意一个方向最大位移量大于2mm或其旋转角度超过2°，剔除该数据；

步骤25：配准，以所述功能磁共振数据头动校正后生成的平均值文件为参考，将所述结构磁共振数据做刚体变换，变换到所述功能磁共振数据的空间里；

步骤26：分割，将变换到功能像的结构磁共振数据分割成灰质、白质、脑脊液，并生成一个将功能空间配准到标准空间的旋转信息矩阵；

步骤27：空间标准化，将所述旋转信息矩阵应用到所述功能磁共振数据里，从而将其自动空间标准化到标准成像模板上，同时标准化后的功能磁共振数据被重采样为3mm各向同性体素；

步骤3：对预处理后的功能磁共振数据进行ppi建模分析，ppi建模分析用于解析感兴趣区bold信号之间的相关性是否受到任务参数的调控，如果是则称感兴趣区间相关性与任务参数存在交互；所述感兴趣区为根据解剖脑图谱模板将整个人脑划分为116个脑区，每个脑区为一个感兴趣区；ppi建模分析则是计算每个感兴趣区的bold信号与任务参数交互的脑区激活信号，采用模型方程如下式：

y＝aβ1+bβ2+abβ3+gβ4+ε

其中，因子a为感兴趣区，因子b为任务参数，aβ1为因子a的主效应项，代表感兴趣区的bold信号，bβ2为因子b的主效应项，代表任务参数相关的blod信号，abβ3为心理生理交互项，即因子a与因子b的交互项，代表受任务参数调控的感兴趣区的bold信号，gβ4为协变量项，ε为误差项，g是头动参数，β1、β2、β3和β4为相应项的参数，y为预处理后的功能磁共振数据；

建立ppi模型的步骤包括：

步骤31：第一次广义线性模型分析，用任务参数b和头动参数g构建的设计矩阵来拟合预处理后的功能磁共振数据，从而构建线性回归方程模型并进行模型估计，估计出任务参数相关的bold信号bβ2；

步骤32：设定任务参数对比条件，同时定义感兴趣区，从限定了对比条件的所述第一次广义线性模型结果中提取感兴趣区的bold信号aβ1；

步骤33：生成所述感兴趣区bold信号aβ1与任务参数b的交互项，将步骤32中提取的所述感兴趣区bold信号aβ1与任务参数b相卷积，生成心理生理交互项abβ3；

步骤34：第二次广义线性模型分析，用心理生理交互项abβ3、所述感兴趣区bold信号aβ1、任务参数b和头动参数g构建的设计矩阵再次拟合预处理后的功能磁共振数据，从而构建线性回归方程模型并进行模型估计，估计出任务参数相关的bold信号bβ2与所述感兴趣区的bold信号aβ1交互后激活的全脑bold信号；

步骤35：步骤31-步骤34为对单个感兴趣区的ppi建模分析过程，对全脑116个感兴趣区的ppi建模分析，则对每个感兴趣区重复步骤31-步骤34，得到116个感兴趣区的ppi分析结果；

步骤4：定义网络节点，采用116个脑区作为网络节点，一个节点对应116个脑区中的一个脑区；

步骤5：定义网络边，网络边代表网络节点之间的连接，网络边的连接值的计算方法为：从步骤34中得到的激活的全脑bold信号中分别提取116个节点对应的信号值，得到一个116×116的功能连接信号矩阵，两两节点之间的信号值则代表这两节点之间的连接；

步骤6：根据所述网络节点和所述网络边，构建基本脑网络；

步骤7：对所述基本脑网络进行图论分析，通过对构建的基本脑网络的网络拓扑属性进行分析从而解析任务状态下的脑网络模式，包括步骤：

步骤71：二值化网络图谱，选定稀疏度阈值，稀疏度定义为一个脑网络中所存在的边的数量与该脑网络可能的边数的最大值的比值，根据稀疏度阈值将所述功能连接信号矩阵二值化；

步骤72：网络拓扑属性计算，包括全局网络测量值和局部网络测量值；

步骤73：网络参数统计分析，对所述全局网络测量值和所述局部网络测量值进行非参数排列检验，以得到与任务参数相关的脑网络属性。

本发明的有益效果在于：

1、本发明针对现有的脑网络构建方法的不足进行了改善，采用心理生理交互(ppi,psychophysiologicalinteraction)方法计算并提取出受任务参数调制的各脑区间bold信号的功能连接。采用该种方法提取出的连接矩阵再用于网络构建，所构建的网络将融合与任务参数交互的信息，为有效反映认知控制与功能磁共振脑信号交互的脑网络研究提供更有效的生理意义。

2、通常情况下，现有的脑网络分析都是直接从预处理后的功能磁共振数据里提取节点的bold信号，并计算节点间bold信号的偏相关系数，该系数值再作为节点间的边的连接值。而本发明的改进之处在于，先对预处理后的功能磁共振数据进行ppi分析，提取受任务参数调制的各脑区间bold信号的功能连接值，再将功能连接值作为脑网络分析中节点间边的连接值，使得所构建的脑网络涵盖了任务参数相关信息。

附图说明

图1是本发明脑网络构建的方法流程图；

图2是本发明功能磁共振数据预处理的流程图；

图3是本发明ppi分析过程的流程图；和

图4是本发明根据ppi分析得到的数据构建脑网络的流程图。

具体实施方式

下面结合附图进行详细说明。

本发明中的刚体变换是指：只有物体的位置(如平移变换)和朝向(如旋转变换)发生改变，而形状不变，得到的变换称为刚体变换。

本发明中的广义线性模型是指：广义线性模型(generallinearmodel,glm)数学表达式：y'＝β'*x'+ε'，其中ε'表示误差。

广义线性分析基于这样一种假设：每个体素上的磁共振数据(用y'来表示，同一任务的时间序列，即bold信号)，是一些未知参数(用x'表示)的线性组合(用β'表示)。这些未知参数是与任务或时间有关但与具体脑区无关的已知参数，它组成的矩阵通常又称为设计矩阵。在经过glm描述之后，分析的对象就发生了转变。原本是对y'的统计分析，现在转向了使用最小二乘法拟合得到的参数β'，再对β'系数做t检验或者f检验，并根据t检验或者f检验统计量和阈值，对β'系数做统计推断。

我们针对现有的脑网络构建方法进行了改善，采用心理生理交互(ppi,psychophysiologicalinteraction)方法计算并提取出受任务参数调制的脑区bold信号。采用该种方法提取出的脑区bold信号所构建的网络将融合与任务参数交互的bold信息，为有效反映任务参数与功能磁共振脑信号交互的脑网络研究提供更有效的生理意义。

图1是本发明脑网络构建的方法流程图。

如图1所示，本发明的一种功能磁共振数据预处理方法，方法包括：

步骤1：采集被试大脑在做任务状态下的磁共振数据，磁共振数据包括结构磁共振数据和功能磁共振数据，其中功能磁共振数据是一个四维的时间序列图。

步骤2：对功能磁共振数据进行预处理，图2是本发明功能磁共振数据预处理的流程图。如图2所示，预处理步骤如下：

步骤21：去除功能磁共振数据前四个时间点的数据。

步骤22：时间层校正，采用辛格插值法对功能磁共振数据进行时间偏移量校正。时间层校正目的是用来校正同一个功能磁共振数据中层与层之间采集时间的差异，使得层与层之间采集时间被校正到同一个时间点。

步骤23：头动校正，将时间层校正后的中间层图作为基准图像，然后通过旋转或平移等刚性变换将时间序列上的其它层图像与第一层图像的位置匹配，然后用内插值算法对磁共振数据进行重新采样。同时生成头动的位置变化和旋转角度变化图，以及头动文件，头动文件包含六列头动参数。

步骤24：减少运动伪影，若每个被试者的功能磁共振数据头动的任意一个方向最大位移量大于2mm或其旋转角度超过2°，剔除该数据。

步骤25：配准，以功能磁共振数据头动校正后生成的平均值文件为参考，将结构磁共振数据做刚体变换，变换到功能磁共振数据的空间里，从而保证功能像与结构像在同一个空间位置。具体效果为将同一个被试的磁共振数据配准到一起，如功能磁共振数据配准到该被试的结构磁共振数据上。

步骤26：分割，将变换到功能像的结构磁共振数据分割成灰质、白质、脑脊液，并生成一个将功能空间配准到标准空间的旋转信息矩阵。

步骤27：空间标准化，将分割后生成的旋转信息矩阵应用到功能磁共振数据里，从而将其自动空间标准化到标准成像模板上，同时标准化后的功能磁共振数据被重采样为3mm各向同性体素。此步的目的是将不同容积及形状的被试的功能磁共振数据放到一个标准空间里，用一个公用的坐标系去描述具体的一个位置。

图3是本发明ppi分析过程的流程图。

如图3所示，对功能磁共振数据进行心理生理交互分析方法，步骤包括：

步骤3：对预处理后的功能磁共振数据进行ppi建模分析，ppi建模分析用于解析感兴趣区bold信号之间的相关性是否受到任务参数的调控，如果是则称感兴趣区间相关性与任务参数存在交互；感兴趣区为根据解剖脑图谱模板将整个人脑划分为116个脑区，每个脑区为一个感兴趣区；ppi建模分析则是计算每个感兴趣区的bold信号与任务参数交互的脑区激活信号，采用模型方程如下式：

y＝aβ1+bβ2+abβ3+gβ4+ε

其中，因子a为感兴趣区，因子b为任务参数，任务参数可以理解为心理参数，aβ1为因子a的主效应项，代表感兴趣区的bold信号，即感兴趣区的神经活性，bβ2为因子b的主效应项，代表任务参数相关的bold信号，即任务参数对应的神经活性，abβ3为心理生理交互项，即因子a与因子b的交互项，代表受任务参数调控的感兴趣区的bold信号，gβ4为协变量项，ε为误差项，g是头动参数，β1、β2、β3和β4为相应项的参数，y为预处理后的功能磁共振数据。

建立ppi模型的步骤为：

步骤31：第一次广义线性模型分析，用任务参数b和头动参数g构建的设计矩阵(即参数x')来拟合预处理后的磁共振数据(即y')，从而构建线性回归方程模型并进行模型估计，估计出任务参数相关的bold信号bβ2。

步骤32：设定任务参数对比条件，同时定义感兴趣区，从限定了对比条件的第一次广义线性模型结果中提取感兴趣区的bold信号aβ1，即因子a的主效应项。

步骤33：生成感兴趣区bold信号aβ1与任务参数向量的交互项，将步骤32中提取的感兴趣区bold信号aβ1与任务参数b相卷积，生成心理生理交互项abβ3。

步骤34：第二次广义线性模型分析，用心理生理交互项、感兴趣区bold信号、任务参数b和头动参数g构建的设计矩阵再次拟合预处理后的功能磁共振数据，从而构建线性回归方程模型并进行模型估计，估计出任务参数相关的bold信号bβ2与感兴趣区的bold信号aβ1交互后激活的全脑bold信号，即全脑激活信号。全脑激活信号表示其与源感兴趣区bold信号之间的功能连接(即相关性)受到任务参数调制或影响。

步骤35：步骤31-步骤34为对单个感兴趣区的ppi建模分析过程，对全脑116个感兴趣区的ppi建模分析，则对每个感兴趣区重复步骤31-步骤34，得到116个感兴趣区的ppi分析结果，及116个感兴趣区的激活信号图。

图4是本发明根据ppi分析得到的数据构建脑网络的流程图。如图4所示，包括以下步骤：

步骤4：定义网络节点，采用解剖脑图谱模板将大脑划分的116个脑区作为脑网络节点，一个节点可以理解为116个脑区中的一个脑区。

步骤5：定义网络边，网络边代表网络节点之间的连接，网络边的连接值的计算方法为：从步骤3中得到的116个感兴趣区的激活信号图中分别提取116个节点对应的信号值，得到一个116×116的功能连接信号矩阵，两两节点之间的信号值则代表这两节点之间的连接。

本发明中节点的选择标准与步骤3中ppi分析中所选取的感兴趣区一致，均为116个脑区aal标准模板，aal全称是automatedanatomicallabeling，all分区是由mni(montrealneurologicalinstitute)机构提供，all模板将整个大脑划分为116个皮质、亚皮质以及小脑感兴趣区，其中，左右大脑各45个、小脑26个。该模板也被用于步骤3的ppi分析，并得到的每个被试的116个脑区与任务参数交互的激活信号图。

步骤6：根据网络节点和网络边，构建基本脑网络。

步骤7：对基本脑网络进行图论分析，通过对基本脑网络的网络拓扑属性进行分析从而解析任务状态下的脑网络属性。

图论分析步骤如下：

步骤71：二值化网络图谱：选定稀疏度阈值，将网络稀疏度阈值s应用到步骤5中生成的功能连接信号矩阵中，稀疏度s定义为一个脑网络中所存在的边的数量与该脑网络可能的边数的最大值的比值。

该方法能够将所有产生的网络标准化，通过应用特定对象的相关系数阈值，使其具有相同的节点数和边数，并且减少组间总体相关强度值的可能偏差效应，从而能够探索相关脑网络组织的组间差异。然而，到现目前为止还没有专门的定义方法来选择一个单阈值，本发明的技术方案不止选择一个单阈值，而是在一个较宽范围的稀疏度水平间反复对每个功能连接矩阵进行阈值化。

稀疏度的阈值范围的选择标准如下：

1)每个阈值化网络的所有节点的平均节点度大于2倍的log(n)(约为4.1289)，其中n＝116，代表着节点的数量。一个节点的节点度定义为连接到该节点的边数。

2)所有被试的阈值网络的小世界属性的标量值σ均大于1.1。小世界网络是一种介于随机网络和规则网络之间的网络，其小世界属性具备条件为：从局部看小世界网络具有高度的集团化特征，而从整体看任意两点间的平均路径长度有能够保持较短。在选定的稀疏度范围内，每个脑网络的最大组成度在110到116之间，最大组成度是指脑网络中所包含最多的节点个数。随后的网络分析按0.01的步长在所选择的稀疏度范围(经验范围值：0.09<s<0.4)反复执行。

步骤72：测量脑网络拓扑属性，包括全局和局部网络测量值。其中，全局网络测量值包括：1)涉及聚类系数cp，特征路径长lp，标准化的聚类系数γ，标准化的特征路径长λ，以及小世界标量σ在内的小世界参数；2)涉及局部效率eloc和全局效率eglob在内的脑网络效率。局部网络测量值包括三个节点中中度测量值：即节点度k，节点效率e，以及中介性b。

其中，聚类系数cp用于量化一个网络的局部互连或聚集程度。聚集程度越高代表局部处理信息的效率越高。特征路径长lp度量一个网络中的任意一对节点之间的距离(即边的数量)，或一个网络的整体连接效率程度。距离越短代表远程效率更高，因为信息通过较少的步骤进行交换。全局效率和局部效率分别衡量一个网络在全局和局部状态进行传输信息的能力。

步骤73：网络参数统计分析，为了得到稳定的任务参数相关脑网络属性，需对所测量到的网络拓扑属性进行非参数排列检验，该分析方法被用于计算每个脑网络矩阵测量值，脑网络矩阵测量值包括小世界属性、脑网络效率和局部中心度测量参数。简而言之，我们首先计算每个脑网络度量的平均值的组间差异。为了检测观测到的组间差异的零假设是否会偶然出现，对于每一个网络参数，我们重新随机分组，并计算随机分组后的平均差异，零假设是指进行统计检验时预先建立的假设。随机过程重复1000次，并采用每个分布的95％作为零假设的单尾t检验的临界值，将一类错误控制在0.05内。值得注意的是，在进行排列检验前，多重线性回归分析用于去除协变量，包括每个脑网络度量的年龄和性别，其中，每个脑网络度量作为独立变量；年龄和性别作依赖变量。

步骤73的作用在于检验步骤72所测量的全局和局部脑网络测量值是否显著，如果显著则可以用于解释该属性与任务相关。

需要注意的是，上述具体实施例是示例性的，本领域技术人员可以在本发明公开内容的启发下想出各种解决方案，而这些解决方案也都属于本发明的公开范围并落入本发明的保护范围之内。本领域技术人员应该明白，本发明说明书及其附图均为说明性而并非构成对权利要求的限制。本发明的保护范围由权利要求及其等同物限定。