功能磁共振成像熵连接工具箱的研制

石莉莉，刘梓君，穆文迪，焦青，曹卫芳，崔栋，郭永新

山东第一医科大学（山东省医学科学院）放射学院，山东泰安271016

前言

目前，基于血氧水平依赖的功能磁共振成像（functional Magnetic Resonance Imaging, fMRI）技术在脑神经活动研究领域起到举足轻重的作用。通过该技术可对大脑整合功能进行研究，分为功能连接和效应连接，前者是指相对距离较远的脑区神经活动之间的关系，可体现不同脑区间的统计相关性，不具方向性。后者是指脑区神经元活动之间产生的有向性影响，具有方向性。常用的效应连接方法有动态因果模型（Dynamic Causal Modeling, DCM）［1］、Granger 因果分析（Granger Causality Analysis, GCA）［2］及熵连接［3］等。熵连接是一种基于信息论的方法，利用改进的传递熵来测量方向函数连通度的强度，该方向函数连通性定义为熵连接，可描述脑区信号间的因果关系。Zhang等［4］利用熵连接对单侧感音神经性耳聋患者在静息态下脑区间的效应连接情况进行观察，发现与正常听力个体相比，左侧长期单侧感音神经性耳聋患者左侧初级听觉皮层脑区表现增强的同步输出熵连通性，左侧感音神经性耳聋患者比右侧感音神经性耳聋患者表现出更显著的熵连接性变化。

本研究基于科研实际需要，在Windows 7.0 系统下，利用MATLAB 平台中图形用户界面设计环境开发了一个熵连接计算工具箱（Entropy Connectivity Analyze Toolbox,ECAT），用于计算脑区间fMRI信号的熵连接。

1 熵连接算法

设有时间序列x和y，xt和yt分别是x和y在t时刻的值。从x上取N个样本，在时刻t处x的增量为Δxt=xt+1‐xt，时刻t+1处y的增量为Δyt+1=yt+2‐yt+1，t= 1, 2,…,t‐2。设n和m分别表示Δxt和Δyt+1之间同步和异步变化的总步长，同步变化指时间序列x和时间序列y均增加（减少），异步变化指时间序列x和时间序列y增量变化方向相反（一个增加同时另一个减少）。若t时刻Δxt和Δyt+1同时大于零或同时小于零，则n加1，由表示同步变化的贝叶斯概率［5］。反之，当Δxt大于零，Δyt+1小于零或Δxt小于零，Δyt+1大于零时，则m加1，由m/(N‐2)表示异步变化的贝叶斯概率。rxy表示x和y之间的皮尔森相关系数，则x与y之间的同步熵连接TSx定义为：

x与y之间的异步熵连接TAx定义为：

由n1和m1分别表示Δyt和Δxt+1之间同步和异步变化的总步长，则Ps(y/x) =n1/(N‐2)为同步变化的贝叶斯概率，Pa(y/x) =m1/(N‐2)为异步变化的贝叶斯概率，则y和x之间同步熵连接TSy定义为：

y和x之间异步熵连接TAy定义为：

则由x指向y的方向同步熵连接TSx→y定义为：

由x指向y的方向异步熵连接TAx→y定义为：

TSx→y≥0 或TAx→y≥0 表示熵连接的方向是由x到y。相反，TSx→y＜0 或TAx→y＜0 表示熵连接的方向是由y到x。远离某个区域的熵连接称为该区域的输出熵连接，指向某个区域的熵连接称为该区域的输入熵连接。

2 功能实现

本研究采用MATLAB 语言编写熵连接计算函数，用于在GUI中进行调用。ECAT 工具箱主界面见图1，有中文（图1a）与英文（图1b）两种版本，由熵连接计算及统计分析两部分组成。熵连接计算模块对用户输入的感兴趣区域（Region of Interest,ROI）的时间序列或Nii 图像文件进行处理，从而得到脑区与脑区之间的熵连接强度与方向。统计分析部分可进行熵连接的单样本t检验和双样本t检验。

2.1 熵连接计算

该模块包括选择输入栏、选择脑区栏与输出栏，用户可以通过两种方式输入数据。

若用户已经生成了脑区ROI时间序列，可以直接点选“选择ROI时间序列”（图2），然后选择已生成的ROI时间序列文件，该时间序列可以来自静息态fMRI数据分析工具箱［6］及静息态fMRI数据分析助手［7］等软件。此时工具箱利用uigetfile函数创建标准的对话框，通过交互式操作获得存放ROI时间序列的文件路径，采用load 函数读取ROI时间序列为mat文件，利用for循环和熵连接函数来设置回调函数，计算熵连接数值。同时为了作图需要，用户要在“选择脑区”栏中，点选“自定义”，在文本框中输入ROI编号与名称，如1:ACC，2:PCC。

图1 工具箱主界面Fig.1 Main interfaces of the toolbox

图2 选择ROI时间序列计算熵连接Fig.2 ROI time series selected for the calculation of entropy connectivity

若用户没有生成ROI时间序列，可以直接点击“选择nii 图像文件”，然后点选“Brodmann 模板”或者“AAL 模板”提取ROI（图3a），同时根据编号参考（图3b）在编辑框内输入脑区编号，例如［1,2,3,4,5,6］。工具箱中拥有含有脑区编号、名称及坐标的TXT 文件，这样在作图时就可以显示用户所选择的脑区名称。在实现时，工具箱调用dir 函数读取文件夹下被试者的Nifti文件，利用SPM12［8］工具箱中spm_vol函数读取Nifti［9］文件内容，利用spm_read_vols 函数读取三维图像及脑区模板数据，以得到ROI时间序列［10］，最后利用for循环和熵连接函数，编写回调函数以计算熵连接数值。

为了将生成的计算结果保存到指定的路径下，首先利用Uigetdir 函数对弹出菜单控件设置，用于选定文件存放路径，调用mkdir 函数在选定的文件存放路径下创建多个文件夹，输出文件夹名是结果文件前缀与原文件夹名的组合，调用Save 函数，将计算结果mat文档（图3c）并同时保存为Excel文档。

2.2 统计分析

该工具箱统计分析功能模块利用t检验对多个被试的全脑熵连接进行数据分析，可进行单样本t检验和双样本t检验。将多个被试者全脑熵连接的mat文件存放在同一个文件夹下（图4a），在进行单样本t检验时，首先调用Dir 函数获得指定文件夹下所有被试的熵连接结果mat 文件，存放在结构体数组中，然后进入for循环，利用Load函数依次读取mat文件，采用ttest 或ttest2 函数进行单（双）样本t检验，将名为OneT（TwoT）的mat 文件和Excel 文件保存在指定的路径下（图4b）。

3 帮助文档

帮助文档提供关于工具箱的使用说明，在此可找到工具箱使用过程中涉及的所有功能，使用者只需点击列表栏中的目录，即可出现关于此功能的详细解释（图5），该功能的实现是采用了switch 和case函数。

4 应用实例

本研究利用该工具箱，对2 型糖尿病（Type 2 Diabetes Mellitus,T2DM）［11］患者及正常被试的fMRI数据进行分析。患者来自山东省泰山疗养院，共15例，其中男7例，女8例，平均年龄（56.60±7.37）岁。健康人被试来自泰安市社区志愿者，共16例，其中男6例，女10例，平均年龄（56.80±9.94）岁。该研究经过山东省泰山疗养院伦理委员会批准，所有研究对象均对此次研究知情同意。采用GE 1.5T 磁共振成像仪进行静息态血氧水平依赖（BOLD）数据的采集，采用DPABI软件进行BOLD信号的预处理［12］。

图3 选择Nii图像文件计算熵连接Fig.3 Choosing Nii file to calculate entropy connectivity

图4 统计分析Fig.4 Statistical analysis

图5 帮助文档Fig.5 Help documents

为了显示不同脑区之间的熵连接，工具箱可生成脑区熵连接的示意图。对于每一个脑区，它与其他脑区之间都呈现出4种熵连接性，分别是同步输出熵连接、异步输出熵连接、同步输入熵连接、异步输入熵连接。当要对两组人的上述4 种熵连接性进行比较分析时，程序首先利用get 函数读取脑区选择模块中编辑框内容，获得与脑区编号相对应的脑区名称与坐标，然后进入for 循环，调用plot 函数，依据双样本t检验结果，显示两组被试间具有显著性差异（P＜0.05）的熵连接性数值，得到脑区间熵连接示意图（图6），同时利用text 函数在示意图中对脑区名称进行说明。熵连接的方向用箭头表示，熵连接强度利用Arrow 函数通过线段粗细加以表示。最后采用saveas函数，将脑区间熵连接示意图保存到与熵连接计算结果相同的文件夹下。

图6 脑区间同步输入与同步输出熵连接差异示意图Fig.6 Entropy connectivity differences of synchronous inputs and outputs between brain regions

本研究中，选取8 个脑区为ROI，分别是内侧前额叶（Medial Prefrontal Cortex,MPFC）（编号7）、后扣带回（Posterior Cingulate Cortex, PCC）（编号35）、双侧海马（Hippocampus）（编号37, 38）、双侧杏仁核（Amygdala）（编号41, 42）、双侧角回（Angular Cortex）（编号65,66），提取这8 个脑区的fMRI信号，计算脑区间的熵连接，经过两组间比较，将两组间具有显著性差异的熵连接进行显示（图6），图中箭头粗细表示熵连接的强度，箭头方向表示熵连接方向。由图6可见，与正常健康人相比，T2DM 患者的脑区间同步以及异步熵连接具有显著差异，比如从左侧海马的同步输入到MPFC的熵连接升高，而杏仁核接收从MPFC传来的同步熵连接降低（图6a），这些表明脑区间神经同步脑波动协同性的升高或降低。左侧角回传到左侧海马和左侧杏仁核的输出熵连接增强，传到MPFC 和右侧角回的输出熵连接降低（图6b），表明脑区间神经活动的相互抑制性增强或者降低。

5 讨论

本研究从科研实际需求出发，在MATLAB 环境下开发了一个可计算时间序列间的熵连接的工具箱，并可实现对输出结果进行单样本及双样本t检验，以及结果的图形显示。

本工具箱中的熵连接代码，可用于计算fMRI信号脑区时间序列之间的效应连接。脑功能连接反映不同脑区间的统计依赖性（相关或者一致性）［13］，但无法说明一个脑区的神经活动如何影响另一个脑区。效应连接可说明不同脑区时间序列间的因果关系，如GCA［14］和DCM［15］都已经得到了广泛的应用［16‐18］。本研究中的熵连接源自于传输熵［19］，不假设任何的特定模型，而是基于假设静息态fMRI信号是随机的，可近似认为是马尔科夫过程。熵连接将两个脑区之间的因果关系细化到4 个方面，即同步输入、同步输出、异步输入及异步输出。脑区间增强的同步熵连接表示脑区间协同关系增强，同步熵连接减弱表示脑区间协同关系降低。异步熵连接表示区域间信息流是异步的波动，增强的异步熵连接脑区间相互抑制增强。本研究中代码运行结果与算法提出者的网络公开代码运行结果一致［20］。

工具箱界面友好，操作简单。用户无需了解过多的操作知识，只需掌握普通的电脑操作，便可实现程序运行操作。工具箱充分考虑到用户需求，输入多样化，输出结果便于用户进一步处理。工具箱可接收用户自己已经生成的ROI时间序列，也可以接收预处理后的图像文件，通过现成模板为用户生成ROI时间序列，从而计算熵连接数值。另外，计算结果可保存为Excel 格式文件，同时保存成MATLAB 下的mat 数据文件，便于用户进行进一步处理。同时，可实时给出熵连接示意图，显示熵连接强度与方向，并标注脑区名称。该结果示意图可保存为MATLAB下的figure文件，便于用户进行进一步编辑。

MATLAB 语言具有很好的开放性与易操作性，便于工具箱的后期完善与更新。在MATLAB 环境下只需拖拽需要的控件到界面的相应位置，即可完成GUI的界面布局。后期若需要扩展功能，只需添加新的计算模板，无需对界面做大规模调整。

综上所述，本研究的工具箱可供科研人员计算fMRI数据但不局限于fMRI数据的效应连接，其输入可以是任何的时间序列或者是一系列数据，因此具有广阔的应用前景，同时该工具箱也可作为教学工具向学生展示fMRI数据处理过程，界面友好，操作方便，具有一定的实用价值。