张笑非，钟宁，杨阳，黄佳进，陈建辉，马莲芳

1.北京工业大学信息学部，北京 100124；2.江苏科技大学计算机学院，江苏镇江 212003；3.前桥工科大学生命科学与信息工程系，日本群马前桥 371-0816；4.北京林业大学心理系，北京 100083； *通讯作者 钟宁 zhong@maebashi.ac.jp

心算是在不借助其他工具和设备帮助的前提下，完全通过大脑功能完成基本的算术过程。既往研究表明，心算会显著激活后顶叶的前中央回和小叶，即处理和工作记忆的中心[1]；甚至不同心算的亚型在大脑区域被激活的程度也存在差异。有研究表明，左侧额下回、背外侧前额叶皮质的中间部分和辅助运动区在减法过程中比加法过程更活跃[2]。心算会导致双侧顶内沟和枕下回的激活更强烈；且与算术概念处理活动相比可更强烈地触发顶叶和枕叶之间的功能连接[3]。fMRI 相关研究表明，在数字与非数字顺序判断中，顶内沟具有相似的激活[4]。顶内沟是表示数值定量语义的关键节点[5]。两位数字的加法和减法均依赖于顶下小叶内区域、前突神经以及枕下回和某些皮质下结构。但加法更加依赖额顶皮层回路；而减法则更依赖顶皮层下回路[6]。脑功能网络分析在精神疾病、高级认知等领域得到广泛应用。多项研究对网络的功能连接变化展开了深入探索[7-10]。目前，对于心算认知状态下脑功能网络以及不同心算状态下脑功能网络变化的研究较为有限，更为深入的研究将有助于理解人脑在进行心算时的组织机制和动力学。本研究拟通过脑图Dosenbach-160 构造受试者在心算试验两种认知 状态下的功能网络，计算相应的网络指标，并在两种心算认知状态间进行对比，从而观察脑功能网络在心算过程中的动态调整。

1 资料与方法

对心算试验中21 组成年人tfMRI 数据进行基于脑图Dosenbach-160 的图论分析。初始假设加法和减法认知状态下脑功能网络的拓扑结构存在较大差异，且网络指标存在显著差异[11]。为验证这一假设，首先为两种心算认知状态构建功能连接网络；然后计算网络拓扑的定量指标，包括全局效率、局部效率、聚类系数、度中心性等；最后对两种心算认知条件下的脑功能网络指标进行统计检验。研究整体流程见图1。

图1 研究的整体流程

1.1 受试者与试验设计 21 名受试者参与此项心算试验。其中男12 名，女9 名，平均（25.76±3.66）岁。受试者均为右利手，具有相同的教育背景，正常或矫正的正常视力，无神经或精神病史。在试验开展前，向每名受试者解释试验过程可能产生的自然反应。所有受试者均签署知情同意书，本研究获得首都医科大学宣武医院伦理委员会批准。

试验中的刺激条件为简单加法和简单减法，即2 位数不产生进位的加法及2 位不包含借位的减法。试验期间，受试者通过观看图片的形式接收刺激；并在每次试验中，依次为受试者呈现4 张图片，其中第1 张和第2 张为本次试验的2 个操作数，第3 张为操作符（加号或减号），第4 张为推荐答案。受试者可通过触摸按钮方式判断第4 张图片中的结果是否正确。试验中要求受试者尽可能快速和准确地点击。记录每次试验的反应时间和正确度。此外，试验还包括一项数配试验和睁眼休息。数配试验要求受试者判断推荐答案是否与某项操作数一致；睁眼休息区块的持续时间与其他类型区块相同，旨在使受试者的认知状态恢复至相对初始化的水平。

由于心算试验的持续时间较长，故整个过程分为2 次试验节段。如图2A 所示，每次试验节段的持续时间为552 s，受试者被允许在2 次试验节段之间进行片刻的休息。具体时间根据实际情况调整。如图2B所示，每次试验节段中包含23 个试验区块，每个区块的持续时间为24 s。加法、减法以及数匹区块在每次试验节段中分别出现4 次，其间会插入休息区块。因此，休息区块仅出现在偶数区块上。因此在加法、减法和数匹区块之后，受试者具有固定的休息时间。类似地，加法、减法和数匹区块仅出现在奇数区块上，并且出现的随机顺序将根据某些规则确保平衡。如图2C 所示，每个区块内进行连续4 次同类型试验。如图2D 所示，每个试验中的操作数显示持续0.25 s，操作符显示持续0.5 s，推荐答案显示持续2 s。操作数、操作符和推荐答案之间的黑色背景持续0.5 s，推荐答案后的黑色背景持续1.5 s。

图2 心算试验范式

1.2 数据采集和预处理 应用3T MAGNETOM Trio Tim（西门子医疗系统）获得12 通道相控阵头部线圈。试验中限制受试者的头部运动并减少扫描期间产生的噪音，并分别使用T1 加权3D 磁化准备快速梯度回波（magnetization prepared rapid acquisition gradient echo，MPRAGE）序列和回波平面成像（echo planar imaging，EPI）序列获取结构影像和功能影像。fMRI 数据采集使用T1 加权3DMPRAGE 序列：TR 1600 ms，TE 3.28 ms，TI 800 ms，FOV 256 mm×256 mm，翻转角9°、体素1 mm×1 mm×1 mm。EPI 序列：TR 2000 ms，TE 31 ms，翻转角90°，FOV 240 mm×240 mm，矩阵64×64。收集功能图像。获得30 个轴向切片，厚度4 mm，间隙0.8 mm。

fMRI 数据的预处理步骤，第1 步：切片定时校正，其中中间切片用作参考切片。第2 步：重新对准或运动校正，旨在解决时间序列中图像位置的不匹配。第3 步：空间归一化，目的是在标准空间坐标中对受试者间进行平均和激活脑区坐标的报告。本研究采用蒙特利尔神经学研究所（Montreal Neurological Institute，MNI）标准空间，并将所有脑影像重新采样为3 mm×3 mm×3 mm。第4 步：平滑，旨在通过去除高频信息增加信噪比，对齐所有结构以进行受试者间的平均，并在使用随机场理论时提高统计的有效性。此步骤中使用8 mm 半峰宽（full width at half maximum，FWHM）各向同性高斯核进行卷积计算。

由于本研究探讨心算时脑功能连接神经机制，故不仅需要考虑血氧水平依赖（blood oxygenation level dependent，BOLD）信号的频带和脑图谱的选择，还需考虑根据试验条件对fMRI 影像进行重新排序。由于tfMRI 的扫描时间为2 s，故在为受试者构建功能网络时，为消除BOLD 信号中的线性漂移，处理上进行频段为0.008～0.5 Hz 的带通滤波。滤波操作利用软件REST 完成（http://restfmri.net/forum /index.php）[12]。

1.3 脑功能网络构建 为从tfMRI 数据相应体素中提取信号脑功能网络，本研究使用图3中描绘的Dosenbach-160 脑图。该脑图的160 个感兴趣脑区（ROI）通过元分析获得，每个ROI 被唯一地分配给Dosenbach-160 脑图的6 个功能子网之一，其默认网络、前额叶-顶叶网络、扣带回-岛盖网络、感知-运动网络、枕叶网络、小脑网络分别包含34、21、32、33、22、18 个ROI。每个ROI 时间序列的提取以其MNI坐标为中心，对该中心5 mm 半径内所有体素提取的BOLD 信号作平均，从而得到每个ROI 在心算过程中的时序数据。由于整个时间序列涉及多个心算认实验中的区块，故需要在时间上将加法和减法认知状态下的时序数据分别进行抽取。

图3 Dosenbach-160 脑图及其功能子网络

由于试验设计导致相同类型的试验条件区块被分散到试验节段的不同阶段。因此，在分析脑功能连接时，需要考虑如何处理相同类型的试验条件在时间上不连续的区块。Cohen 等[13]提出了2 种不同的区块组合方法，即区块级功能连接和拼接式功能连接。第1 种方法是单独计算每个区块内的时间序列相关性，然后对同一类型试验条件的所有区块的相关矩阵求平均，以得到每个条件的单个平均时间序列相关性；第2 种方法是对每个区块的BOLD 信号进行标准化，以最小化由于信号强度的变化或跨区块的可变性而引起的伪相关性的可能，然后将同试验认知状态类型的不连续区块的BOLD 信号进行拼接。鉴于fMRI 数据在预处理和过滤期间已经标准化，且每个单独区块中再次归一化的效果未知，故本研究在不对每个区块的BOLD 信号作区块内标准化的前提下，采用拼接式功能连接方法计算两种心算认知状态下的相关矩阵。

1.4 脑功能网络的二值化 为利用图论中的拓扑指标对心算时脑功能网络进行度量，需要将Pearson 相关系数矩阵转换为二值矩阵。由于Pearson 相关系数矩阵所有的值构成了一个数值区间，因此可使用阈值确定对应于矩阵中的每个单元的边是否存在。此外，考虑到边为负值的内涵难以解释，因此在后续计算过程中值为负数的边一律认为不存在。Pearson 相关系数矩阵的边若大于或等于当前阈值时表示边存在，二值矩阵相应的单元为1，反之则为0，如此获得的二值矩阵上可进行多种网络拓扑指标的计算。为将多名受试者用于之后的统计检验，因此阈值的选择是根据图中边的密度计算，及整个脑功能网络存在的边占最大可能边数的百分比。本文使用每名受试者所有正值边构成的图密度（1%～100%）和1%的步长获得功能脑网络的二值矩阵。

图4显示一名受试者在完成减法试验的过程中脑功能网络对应的Pearson 相关系数矩阵和连接组。右侧是在正值边构成的图密度（1%～100%）上取阈值后产生的相应二值矩阵。由于研究仅参考正值边，故当图密度为100%时，二值矩阵中不为1 的边即为负值边；后续的网络指标计算则会在所有图密对应的二值矩阵上展开计算。

图4 脑功能网络的二值化过程

1.5 图论分析计算 多种基于图论的度量可用于研究脑功能网络的拓扑性质。本研究采用网络效率（全局效率、局部效率）[14]、聚类系数[15]和度中心性[16]。假设K 是脑功能网络图密度的阈值集合，tk∈T，i=1，2，3，…100，|T|=100，当取ki作为阈值时，脑功能网络对应的Pearson 矩阵转换得到二值矩阵Mk。若k=10，所得二值矩阵中值排在最高10%的正值边保留，其余的正值边与非正值边均认定不存在。但是所有的点均存在，即使该点没有任何边与之连接。这样得到的图这里称为G=（V，E），V 是所有脑区构成的节点集合，E 是所有存在的边的集合，其包含的节点数量为N=|V|。由于这里使用了Dosenbach-160 脑图且接受度为0 的节点，所以在不同的图密度阈值下的N 均为160。下面对图论分析中用到的4 项指标进行说明。首先是全局效率，它是一种可以评价整个网络信息交换率的衡量标准，其定义如公式（1）所示。

其中Di,j表示在图G中从节点vi到vj的最短路径距离。

节点的局部效率就是由该节点的邻居节点所构成网络的全局效率，反映出当前节点的邻居网络的通信效率。而整个网络的局部效率是对所有节点的局部效率进行平均，反映了整体网络的邻居网络通信效率水平。节点vi∈V 的邻居节点构成的集合Mi，对应的图g的全局效率定义如公式（2）所示。

其中di,j表示在图g中从节点vi到vj的最短路径距离。整个图G的局部效率是所有节点邻居网络的全局效率值的均值，定义如公式（3）所示。

网络的聚类系数和特征路径长度常用于量化小世界属性，本文主要采用聚类系数。同样对于特定图密度下产生的二值矩阵所对应的图G，节点vi的局部聚类系数定义如公式（4）所示。

同样，Mi是节点vi的所有相邻节点构成的集合，ejk是存在于节点vj和节点vk之间的边，vj和vk是节点vi的两个不同相邻节点。分子是存在于节点vi相邻节点之间边的数量，|Mi|是节点vi相邻节点的数量，分母是节点vi相邻节点之间可能存在边的最大数量。整个图G的聚类系数是对所有节点的聚类系数求平均得到，如公式（5）所示。

其中Ci是节点vi的局部聚类系数。

度中心性是图论中的一个常见指标，可用来计算给定脑区与其他脑区存在的连接数量。节点vi∈V 的 局部聚类系数定义如公式（6）所示。

同样，Mi是节点vi的所有相邻节点构成的集合，|Mi|是节点vi相邻节点的数量。整个图G的度中心性是对所有节点的度中心性求平均得到，如公式（7）所示。

2 结果

如图5所示，每个子图是两种心算认知状态下脑 功能网络在网络稀疏度定义域上21 名受试者图论度量指标的均值（折线）和标准差（阴影），以及两种心算认知状态在每个网络稀疏度上存在的差异性（红色叉）。图5A、5C、5E、5G 分别为对比两种心算认知状态下Dosenbach-160 脑图全局网络的全局效率、局部效率、聚类系数、度中心性的结果；图5B、5D、5F、5H 分别为对比两种心算认知状态下前额叶-顶叶网络的全局效率、局部效率、聚类系数、度中心性的结果。

图5 两种心算认知状态间统计对比结果

图5A 显示，两种心算认知状态下Dosenbach-160脑图全网的全局效率仅在2 个稀疏度上存在显著差异，由于共考察了100 个网络稀疏度，因此显著性差异占比为2%。图5C、5E、5G 中显示两种心算认知状态下Dosenbach-160 脑图全网的局部效率、聚类系数、度中心性在所有的网络稀疏度上均不存在显著差异，即显著性差异占比均为0%。因此，加法认知状态和减法认知状态在Dosenbach-160 脑区全局网络图论指标上无显著差异。

相应地子Dosenbach-160 脑图的功能子网络、即前额叶-顶叶网络上观察加法和减法认知状态的网络度量是否存在显著性差异。如图5B 和5H 所示，两种心算认知状态下前额叶-顶叶网络的全局效率和度中心性在所有的网络稀疏度上均存在显著差异，即显著性差异占比100%。如图5D 和5F 所示，两种心算认知状态下前额叶-顶叶网络的局部效率和聚类系数存在显著差异的网络稀疏度个数分别为71 个和78 个，即显著性差异占比分别为71%和78%。因此，加法认知状态和减法认知状态在前额叶-顶叶网络图论指标上存在显著差异。

3 讨论

本研究发现，人脑在进行心算认知过程中，脑功能网络的拓扑属性在全局并未发生显著改变，但局部的前额叶-顶叶并网络的拓扑属性存在显著变化，表明在心算认知过程中，前额叶-顶叶网络在不断地进行调整、调度、通信以完成或参与完成相应的心算认知功能。成人的算术能力取决于前额叶-顶叶控制网络，其中顶叶区域包括上顶叶和下顶叶，以及包括下额叶回的额叶区域、中额回和左额上回[17-20]。此外，需要更加深入研究的问题是，加法认知状态下前额叶-顶叶网络的4 项指标均显著高于减法认知状态下对应项目。首先，前者的度中心性显著高于后者，体现了网络节点间边的数量，表明在加法认知状态下，前额叶-顶叶网络内脑区间的互动密度高于减法认知状态。其次，由于加法认知状态下前额叶-顶叶网络内边的数量显著高于减法认知状态，节点间具有更高的通信效率，因此全局效率更高。最后，加法认知状态下前额叶-顶叶网络的聚类系数显著高于减法认知状态，说明前者网络中所有节点的邻居网络内关联密度总体高于后者，从而使得所有节点的邻居网络内通信 效率总体高于后者，故局部效率更高。

针对本研究提出的假设，人脑在进行心算认知活动时，加法和减法认知状态下脑功能网络的拓扑结构存在较大差异，且网络指标存在显著差异。通过图论分析和统计分析，发现该假设在采用的脑图Dosenbach-160 的全网上并不成立；但进而发现假设在Dosenbach-160 定义的子网前额叶-顶叶网络上成立，这一结果与既往多项研究结果一致[21-23]。Kesler 等[21]研究发现，受试者在自适应数感训练后数学技能水平提高的同时，前额叶-顶叶的活动也发生了变化，其中双侧顶叶激活明显增加，额叶纹状体和中间颞叶激活降低。Liang 等[22]研究发现，数字归纳推理属于一般化的心算过程，前额叶-顶叶网络在数值归纳推理中更活跃，在规则识别过程中需要中间表示和长期陈述性知识之间进行更多的交换。Tsang 等[23]研究发现，前额叶-顶叶的白质弥散属性可有效预测儿童的心算能力。

最后，本研究通过脑功能网络分析研究人脑的心算功能，更通用性地发现需要通过大量样本复现当前结果。当结果确认后，这些发现可以成为纵向研究的基础，从而测试以下问题：①在心算的加法和减法过程中，Dosenbach-160 脑图前额叶-顶叶网络中的21 个脑区是否存在拓扑属性的显著变化；②其他诸如Power-264[24]、Willard-499[25]等脑图中均存在对前额叶-顶叶网络的定义，如通过这些脑图进行同样的分析，是否可获得相一致的结果；③如将心算过程在时间轴上进行窗口化，是否可通过动态网路分析方法发现前额叶-顶叶网络同样存在活跃的变化。