0 引言

人脑是自然界最复杂的系统之一，负责人体的感知觉、语言、思维、情感、运动等各种功能和活动，科研工作者一直致力于使用各种新技术来研究和探索人脑的工作原理和运行机制。近年来基于核磁共振成像(Magnetic Resonance Imaging,MRI)的人脑结构重建技术^[1]因其无创性的检测手段和高精度的人脑建模方法为人们研究脑科学提供了有利的工具。2005年美国印第安纳大学的神经学家Sporns等^[2]提出了基于核磁共振成像的人脑连接组(human connectome)概念及技术，提出了利用图论和复杂网络的理论去研究人脑网络，尤其是大尺度(marcoscale-大脑脑区)网络^[3]的结构与功能特征。大脑在处理一个具体的行为时，是由多个紧密联系的脑区群体共同完成的，而这些脑区群体的存在则说明脑网络中具有模块结构^[4]，即模块内的节点之间联系紧密，而模块间的连接较少。

目前复杂网络模块化算法包括Kenighan-Lin算法^[5]、谱平分法^[6-7]、GN(Girvan Newman)算法以及Fast Newman算法等，这些算法都是以二值网络为对象的。Girvan等^[8]于2001年提出基于删除边的分裂式模块化算法；Newman^[9]于2004年提出基于凝聚思想的快速算法；Clauset等^[10]在Fast Newman算法的基础上，利用贪心算法思想进行改进。在上述研究基础上，很多学者应用二值模块化算法对人脑网络展开了研究，Meunier等^[11]研究了年轻和年老的静息态脑功能网络的模块结构，结果表明年轻和年老的人脑均具有模块结构，但是有所不同；He等^[12]研究了健康年轻的静息态脑功能网络的模块结构，得到运动体感听觉区、视觉处理、注意力处理、默认网络以及子皮质系统这五个模块。Hagmann等^[13]使用扩散光谱成像获得了大尺度上的人脑结构网络，利用谱平分法证实了在人脑皮质网络中存在模块化结构，划分出六个模块；van den Heuvel等^[14]使用带模块度的谱平分法对人脑结构网络进行模块划分，得出左前脑、右前脑、左后脑、右后脑四个对称的模块。

但是二值人脑网络是以脑区之间的结构连接与否构建的，以此为基础来辨识模块并不足以完整地体现出人脑结构的特点。因为脑区中包含的白质纤维束连接、皮层厚度等生物指标对脑区模块的形成具有重要的影响，例如白质纤维束连接负责协调脑区之间的运作，是大脑皮层灰质脑区间信息传递的物质基础，而结合这些生理指标所构建的加权人脑结构网络能够体现大脑的更多特征信息，因此对加权人脑结构网络的模块结构进行辨识将更有利于帮助人们完整的理解大脑结构和功能特点。

1 基于MRI的加权人脑结构网络重构 1.1 数据采集及人脑结构网络重构

本文共采集了60个健康人的脑部MRI数据，由美国通用电气公司制造的1.5Tesla 双梯度磁共振成像仪完成，使用头部正交线圈。弥散张量成像(Diffusion Tensor Imaging,DTI)数据采用单次成像平面回波成像序列，共施加15个非共线弥散敏感梯度。DTI主要扫描参数如下：TR=11000ms，TE=74.7ms，采集矩阵为128×128，视野为240mm×240mm。b值为1000s/mm²；同时采集一组b0数据，b=0s/mm²。每个方向上冠状面采集35个层面，覆盖全脑。层厚为4mm，层间距为0.0mm。

本文采用瑞士洛桑理工学院开发的开源脑软件平台Connectome Mapper pipeline^[15]实现了人脑结构网络的重构，其数据处理和重构过程如下：

1) 采集MRI(T1) 和DTI数据；

2) 对T1与DTI原片进行误差校准，以消除采集过程中的机械误差、头动误差等；

3) 对原片数据进行格式转换并进行数据预处理，包括涡流矫正、平滑、去除非脑组织部分等；

4) 使用Freesurfer模板对大脑皮层进行分割，本文采用83个脑区的分割方案^[16]；

5) 利用白质纤维追踪术(tractography)^[17]对白质纤维进行追踪和重建，同时通过纤维过滤去除杂质和过长过短的纤维组织，得到全脑白质的纤维跟踪结构；

6) 将每一个脑区定义为人脑结构网络(下文简称脑网络)中的一个节点，两两脑区之间的白质纤维束连接作为脑区之间相互联系的度量，并将其定义为脑网络中节点间的边。

经过上述处理过程所得到人脑的白质纤维束连接信息和数据，采用图论知识，可以对其作如下具体表述：

1) 每个数据都由2个集合构成：节点集合Nodes和边集合Edges。

2) 节点集合为Nodes={1,2,…,83}，每个节点即代表每个脑区。

3) 边集合为两两脑区之间的白质纤维束连接数目的集合Edges={M₁₁,M₁₂,…,M_ij,…,M_nn}，n=83，i=1,2,…,83，j=1,2,…,83，M_ij代表脑区i和脑区j之间的白质纤维束连接数目。

4) 如果脑区间无白质纤维束连接，则M_ij为0；如果脑区间存在白质纤维束连接，则M_ij大于0。脑区之间的白质纤维束连接数据可用矩阵M来描述：

$M=\left[ \begin{matrix} {{M}_{11}}&{{M}_{12}}&\cdots &{{M}_{1n}} \\ {{M}_{21}}&{{M}_{22}}&\cdots &{{M}_{2n}} \\ \vdots &\vdots &{}&\vdots \\ {{M}_{n1}}&{{M}_{n2}}&\cdots &{{M}_{nn}} \\ \end{matrix} \right]$

(1)

1.2 组平均网络的构建

在本文研究中，采集并重构了60个健康被试的人脑网络，按照人脑网络研究的常用方法，为了得到具有概括性、普适性的白质纤维束连接数据，将这60个数据看成一组，并作组平均处理^[14]，得到了本文实验中所使用的二值及加权脑网络数据。组平均处理思想是在不丢失重要信息的情况下，保证组平均数据里两两脑区间的白质纤维束连接的存在性和有效性。在这里引入有效存在率的概念，即组平均数据中任一连接(不包含脑区自连接)存在的充分必要条件是该组中至少有75%的个体数据都拥有这一连接^[14]，在本文实验中是指至少在45个个体中存在的连接。具体处理步骤如下：

1) 首先去除脑区的自连接，60个数据各自对应的矩阵M里对角线元素置为0。

2) 检查每对脑区间白质纤维束连接的有效存在率是否大于75%。

3) 仅保留60个矩阵M里那些有效存在率大于75%的白质纤维束连接的数量值，剔除其余的纤维束数量值。

4) 计算组平均数据，每对存在的连接，其数量值根据它对应的数据里的连接数量值来计算平均值。

组平均网络为下一步进行网络的阈值化提供了可靠的数据基础。

1.3 加权脑网络的阈值化

为了进行模块结构的研究，需对脑区间的白质纤维束连接的数目设定阈值，从而去构建二值矩阵和加权矩阵。阈值化的具体步骤如下：

1) 假定阈值范围为[T₁,T₂]，如果矩阵M中的元素满足T₁≤M_ij≤T₂，那么将M_ij的值置为1，表征脑区间的连接关系，反之则置为0，表示脑区间无连接。则矩阵M可表示为：

$M=\left[ \begin{matrix} 0&{{M}_{12}}&\cdots &{{M}_{1n}} \\ {{M}_{21}}&0&\cdots &{{M}_{2n}} \\ \vdots &\vdots &{}&\vdots \\ {{M}_{n1}}&{{M}_{n2}}&\cdots &0 \\ \end{matrix} \right]$

(2)

其中${{M}_{ij}}=\left\{ \begin{array}{*{35}{l}} 1,{{T}_{1}}\le {{M}_{ij}}\le {{T}_{2}} \\ 0,{{M}_{ij}}<{{T}_{1}}\ or\ {{M}_{ij}}>{{T}_{2}} \\ \end{array} \right.$。

2) 阈值T依据网络的连通性和小世界性^[18]来选择。连通性原则是保证网络中所有的节点都有边连接，满足式(3) :

$links>(n-1)/2$

(3)

其中：links为网络中的总边数；n为网络总节点数。

本文实验所构建的脑网络中节点总数n为83，由式(3) 可知，网络的总边数links必须大于41。

小世界性原则是因为脑网络是复杂网络，它相对于随机网络具有较大的集群系数和近似的特征路径长度，满足以式(4) ：

$\left\{ \begin{array}{*{35}{l}} \gamma ={{C}_{net}}/{{C}_{random}}\gg 1 \\ \lambda ={{L}_{net}}/{{L}_{random}}\approx 1 \\ \end{array} \right.$

(4)

其中：C为网络聚类系数；L为网络平均特征路径长度；脚标net表示脑网络；random表示随机网络。

3) 总的阈值选择依据是保证脑网络连通性和小世界性的同时去掉更多的伪连接。这里先假定阈值T₁为1，T₂为脑区间白质纤维束连接数目的最大值，阈值调整的增减量为1。逐步增大T₁，减小T₂，同时检测网络的连通性和小世界属性是否还能得到保证。经过实验，对本文所使用的人脑网络数据，当T₁高于28，或T₂低于212时，脑网络的完整性和小世界属性无法得到保证，所以本文最后选定的阈值范围为[28,212]。

1.4 二值脑网络和加权脑网络的构建

二值脑网络中各个脑区之间的连接关系用矩阵B来描述：

$B=\left[ \begin{matrix} 0&{{B}_{12}}&\cdots &{{B}_{1n}} \\ {{B}_{21}}&0&\cdots &{{B}_{2n}} \\ \vdots &\vdots &{}&\vdots \\ {{B}_{n1}}&{{B}_{n2}}&\cdots &0 \\ \end{matrix} \right]$

(5)

其中${{B}_{ij}}=\left\{ \begin{array}{*{35}{l}} 1,28\le {{M}_{ij}}\le 212 \\ 0,{{M}_{ij}}<28\ or\ {{M}_{ij}}>212 \\ \end{array} \right.$。 加权脑网络中各个脑区之间的白质纤维束连接权重关系可用矩阵W来描述：

$W=\left[ \begin{matrix} 0&{{W}_{12}}&\cdots &{{W}_{1n}} \\ {{W}_{21}}&0&\cdots &{{W}_{2n}} \\ \vdots &\vdots &{}&\vdots \\ {{W}_{n1}}&{{W}_{n2}}&\cdots &0 \\ \end{matrix} \right]$

(6)

其中${{W}_{ij}}=\left\{ \begin{array}{*{35}{l}} {{M}_{ij}},28\le {{M}_{ij}}\le 212 \\ 0,{{M}_{ij}}<28\ or\ {{M}_{ij}}>212 \\ \end{array} \right.$。

用色块图表示矩阵B和W，可得到如图 1所示的脑网络矩阵图。图 1的横纵坐标均为83个脑区的标号，其对应的英文名称，如表 1所示。图 1(a)的黑色表示连接系数为1，白色表示连接系数为0；图 1(b)中颜色由浅入深，表示白质纤维束连接数目的不断增大。

2 基于Fast Newman算法的加权人脑结构网络 2.1 应用Fast Newman算法的二值脑网络模块划分

Fast Newman算法是基于贪婪思想的凝聚算法，其中心思想是依次合并有边相连的节点对i和j，并计算每次合并后模块度增量，选取使得模块度值朝着增长最大或者减少最小的方向的节点对进行合并，直到整个网络合并为一个模块。二值网络的模块度计算公式为$Q=\sum\limits_{u\in m}{[{{e}_{uu}}-{{(\sum\limits_{v\in m}{{{e}_{uv}}})}^{2}}]}$，模块度增量的计算公式为$\Delta Q=2({{e}_{uv}}-\sum\limits_{v\in m}{{{e}_{uv}}}*\sum\limits_{u\in m}{{{e}_{uv}}})$，这是在假定网络被划分成m个模块的基础上，u和v是其中某两个模块，e_uu表示模块u中的边数占网络总边数的比例，e_uv表示模块u和模块v间的连边数占网络总边数的比例。

将其应用在1.4节的二值脑网络中，其初始化矩阵仅考虑节点之间的有无连接，可得到该网络的模块化树状图如图 2(a)，每一次合并节点后模块度值的变化如图 2(b)。图 2的纵坐标为划分过程中依次合并节点的步数，本文选取上方横虚线下的网络结构为二值脑网络的模块结构，对应着最大的Q值0.4209，此时网络中有4个模块，由图中的竖虚线加以隔开。图 2(a)的横坐标依次为:68,82,66,18,33,42,43,47,51,75,78,73,74,63,64,76,71,70,81,65,69,62,56,67,58, 72,1,4,36,44,8,45,12,49,53,54,39,7,10,3,46,48,80,5,9,83,35,13,14,20,21,15, 61,52,55,11,77,59,38,79,50,60,2,32,23,30,17,31,19,6,34,37,28,29,41,25,26,2 7,40,16,22,24。

2.2 加权脑网络的模块化算法及实现

加权脑网络包含了更丰富的人脑结构信息，由1.4节的加权脑网络可知，边的权重大小表示脑区之间白质纤维束连接数目的多少，数目越大表示脑区之间连接强度更大，联系就更为紧密，反之亦然。而点的权重大小则指所有与此节点相连边的边权总和，权重越大表示脑区的信息越可能传递出去，也表现出该脑区在全脑网络中地位越活跃。因此如果忽略权重而仅仅去分析二值脑网络，那么将会丢失大量包含在权重中的重要信息。

本文以Fast Newman算法的凝聚思想为基础，以1.4节的加权脑网络为对象，对加权脑网络的模块结构划分算法展开研究。算法的基本思想如下：将单个脑区节点与其他脑区节点的全部白质纤维束连接总数作为每个脑区节点的权重，将全脑所有脑区节点两两之间的白质纤维束连接总数作为网络总权重，以任意两个脑区节点之间的权重与网络总权重之间的关系为基础分析构建加权脑网络的模块度计算公式，将模块度增量作为依据决定脑区节点的合并。

为描述算法方便，给出如下定义：

定义1   加权脑网络中每个节点的权重定义为K_i^w：

$K_{i}^{w}=\sum\limits_{j=1}^{n}{{{W}_{ij}}}$

(7)

定义2   网络总权重，即全脑的白质纤维束连接数的总数目定义为L^w：

${{L}^{w}}=\frac{1}{2}\sum\limits_{i=1}^{n}{K_{i}^{w}}$

(8)

定义3   加权脑网络的模块度定义为Q^w：

${{Q}^{w}}=\frac{1}{{{L}^{w}}}\sum\limits_{i,j\in n}{({{W}_{ij}}-\frac{K_{i}^{w}K_{j}^{w}}{{{L}^{w}}})}*\delta ({{m}_{i}},{{m}_{j}})$

(9)

其中：δ(m_i,m_j)表示如果节点i和j属于一个模块，δ(m_i,m_j)=1，反之为0。

定义4   加权脑网络的模块度增量定义为ΔQ^w：

$\Delta {{Q}^{w}}=\frac{2}{{{L}^{w}}}({{W}_{ij}}-\frac{K_{i}^{w}K_{j}^{w}}{{{L}^{w}}})$

(10)

加权脑网络模块化算法实现的具体步骤如下所示：

1) 初始化网络，令加权脑网络中每个节点为一个独立模块，即网络中有83个模块，此时Q^w为0，因为节点i和j不属于一个模块。

2) 计算模块权重K_i^w、网络总权重值L^w。

3) 依次合并有边相连的模块对，并计算合并后模块度增量ΔQ^w。

4) 选取使模块度增量为最大值的模块对进行合并，将这两个模块视为一个模块，则网络中模块数减一。

5) 合并后更新网络，将与i和j相关的行和列相加，按照Q^w的计算方法来计算此时网络的模块度。

6) 重复第2) 步至第5) 步，直到整个网络合并为一个模块。

7) 选取最大模块度值所对应的网络结构。

整个算法完成后，可以得到一个树状图，如图 3(a)，每一步得到的模块度值如图 3(b)。图 3的纵坐标表示算法实现过程中合并节点的步数，那么选择在不同的步数位置断开则可以得到不同的网络结构。理论上来说，模块度的取值范围是[0,1]，而在实际的网络中，模块度值通常在0.3到0.7之间。选择局部最大的模块度所对应的网络结构，即为较好的模块结构^[9]。

图 3中最上方的横虚线下的网络结构对应较好的模块结构，此时有着最大的Q^w值0.5387，网络中有5个模块，由竖虚线加以隔开。图 3(a)的横坐标依次为:17,31,30,32,23,19,2,6,37,16,40,59,25,34,28,18,29,24,22,41,56,33,26,27,11,5 2,55,57,13,20,14,21,62,51,1,45,3,8,12,7,44,39,4,49,36,10,9,35,53,5,77,54,8 3,38,80,15,61,79,46,48,68,82,42,47,63,70,71,64,78,76,73,75,65,69,66,72,50, 60,58,43,74,67,81。

3 实验结果分析与对比 3.1 模块结构组成对比

从模块结构本身的意义来讲，它可以帮助人们发现网络中结构与功能间的潜在联系。由图 2可知，二值脑网络包含4个模块，其中模块一包含了左半脑的额叶下方和颞叶大部分脑区，模块二包含了左右半脑的额叶上部以及其周围小部分脑区，模块三包含了左右楔叶和扣带回等脑区，模块四包含了右半脑的额叶下方和颞叶周围脑区。由图 3可知，加权脑网络包含5个模块，模块一里的缘上回负责精细的协调功能，模块二的海马、岛叶以及杏仁核等属于边缘系统，模块三的中央沟回与人的躯体运动有关系，模块四的额上回、额中回以及丘脑主导着人的联合运动与姿势动作，模块五中包含的距状旁回和枕叶外侧对人的视觉形成有着不可或缺的作用。可以看到，使用本文算法的加权脑网络模块划分结果更符合已知的人脑的生理特征。

3.2 模块结构图分析

按照图 3的加权脑网络的模块结构，将此时的网络结构用每个模块内节点的顺序重新排列，获得新的模块结构矩阵，并以彩色图的形式表示出来，如图 4(a)所示。从图 4(a)中可以明显观察到对角线周围的连接数目多，且颜色较深，表明模块内相互连接数目较密集且强度较大，而对角线以外区域的连接很少且颜色很浅，说明模块间的连接数较稀疏且较弱。对比图 4(b)中所示的二值脑网络模块结构矩阵，可以看出其模块间的连接数目仍然较多，且由于采用二值(有无连接)表示，并不能表现出模块内部连接重要性与否。

将加权脑网络的模块结构以节点连接的形式表现出来如图 5所示，五个不同的圆圈分别代表五个模块的组成节点，节点间连接的强弱用线的粗细程度表示出来，线越粗说明节点间的连接强度越大，而线越细则表明节点间的连接越弱。从图 5也可以看出，每个模块内节点的相互连接的强度很大且相对密集，而模块间的连接强度很弱且较少。

3.3 模块度比较

由模块度的定义可以了解到模块度值越大说明网络的模块化程度较高，网络具有明显的模块结构，故模块度是一个有力的模块评价指标。对于二值脑网络的模块结构来说，较好的模块结构所对应的模块度为0.4209，而加权脑结构网络的较好模块结构对应的模块度为0.5387，故加权脑网络的模块度比二值脑网络的提高了28%。而且加权脑网络所识别的人脑模块结构更加合理，模块内连接强度也更大。如果要将二值脑网络的模块结构也同样划分成五个，可以在图 2中将横虚线向下挪一步，但此时对应的模块度为0.4136，比加权脑网络的模块度又有所降低，所以加权脑网络的模块结构更明显紧凑。

3.4 与现有加权网络模块化算法的比较

对本文数据分别使用Blondel等^[19]提出的探测复杂网络中层次化模块的算法，以及Newman ^[20]提出将模块度与谱向量结合的算法来进行加权脑网络的划分，与本文给出的加权Fast Newman算法在模块数量和模块度值上同时进行比较的结果如表 2所示。从最终的模块划分结果可以看出，尽管获得了相近的模块度，但Newman2006算法是基于谱向量的二分思想，只能划分偶数个模块，因此对脑网络模块结构的识别具有较大局限性且实际的结果与已知人脑的生理模块结构相差较大，而采用Blondel算法划分的模块结构过于零散不够紧凑，所得模块化结构未能体现出模块内部与模块间的特征区分。实验结果对比显示，相对于上述两种加权网络模块化算法，本文算法对加权人脑结构网络的模块化辨识更加适用，并且能够保证模块度的小幅提高。

4 结语

在现有的二值复杂网络Fast Newman模块化算法研究的基础上，本文给出了一种包含权重因素的加权Fast Newman模块化算法，并将其应用于以白质纤维束连接数作为权值的加权人脑结构网络中。实验结果表明，与现有的二值Fast Newman算法相比较，该算法能够更准确合理辨识人脑网络的模块结构。与现有的加权模块化算法相比，本文算法能够获得较高的模块度，也能更准确地反映出人脑的生理结构特征。

尽管相对于二值模块化算法，本文所提加权算法表现出一定的优势，但由于Fast Newman算法所使用的凝聚节点思想是以节点的两两合并为基础展开的，而人脑网络权值的个体差异性较大，权值的细微差别所引起的节点凝聚差别将对后面的节点凝聚产生较大影响，因此下一步的研究重点是如何降低该算法对个体权值差异的敏感性。