基于KL 散度的轻度认知功能障碍患者的个体脑结构网络研究

吕 旻，丁 皓，张雪君，张 敬

（1.天津医科大学医学影像学院，天津300070；2.天津医科大学总医院医学影像学科，天津300052）

阿尔茨海默症（Alzheimer’s disease,AD）是一种 以记忆和认知功能逐渐降低为主要临床表现的中枢神经系统退行性疾病，随着病情发展，脑组织出现不可逆的病理损害[1]。轻度认知功能障碍（mild cognitive impairment,MCI）是介于正常老龄化和阿尔兹海默症之间的疾病状态，其临床表现为认知功能或记忆功能的轻度障碍。影像学表现为海马和内嗅皮层等脑组织出现明显萎缩[2-3]。MCI 患者的临床症状具有不确定性，相当多的MCI 患者病情会出现恶化，进而进入AD 阶段。因此，从全脑水平认识易转化为AD 的MCI 患者的脑结构神经机制迫在眉睫。

近年来，越来越多的脑结构影像研究通过利用脑皮层形态特征（如体积、厚度等）来探明脑皮层间的联系[4]，该方法通过计算脑区间皮层形态的共变程度来衡量脑区之间的联系，进而构建脑结构网络。而且，该共变关系可作为生物学标记去探测AD[5]、精神分裂症[6]和癫痫[7]等神经系统疾病引发的脑皮层异常改变。此方法通过采集大量受试者数据，提取多名受试者脑结构像中的每个脑区形态特征，按照预先划定的脑图谱模板依次计算每个脑区与其它脑区形态特征在跨被试间的Pearson 相关系数，并构建共变关系矩阵，进而反映在组水平的脑区间关系[4]。这种方法仅考虑受试者间群体效应，而忽略受试者间脑结构的个体差异。因此，本文利用基于KL 散度（Kullback-Leiblerdivergence）并结合脑皮层形态特征（灰质体积）的方法构建个体化脑结构网络，并通过图论的方法探究正常老龄人群和易转化AD 的MCI 患者脑结构网络拓扑属性指标的差异。

1 资料与方法

1.1 资料选择 本文所使用的实验数据均来自阿尔茨海默神经影像倡议 (Alzheimer’s disease Neuroimaging Initiative, ADNI) 公开数据库(http://adni.loni.ucla.edu/ )，共纳入348 例受试者MRI 数据，包括正常老龄人群（NC）178 例和易转化为AD的MCI 患者（cMCI）170 例。其中，易转为AD 的MCI 患者组是以随访9 年时间内是否由MCI 转化为AD 作为标准纳入数据。

1.2 MR 成像设备及参数 全部脑结构像MRI 数据均由磁场强度为1.5T 磁共振扫描仪依照ADNI协议标准进行采集[8]，通过MPRAGE 序列扫描获得，并且均经过梯度不均匀矫正等数据预处理，具体过程详见ADNI 官网。扫描前，每位受试者均按照要求认真阅读核共振扫描安全须知并填写检查知情同意书。受试者接受扫描前，需由工作人员告知每位受试者扫描试验的内容、扫描试验所需时间以及其他注意事项等。由于扫描过程中噪音较大，受试者需要佩戴耳塞以减少噪音对自身的干扰作用。受试者仰卧于扫描床上，扫描人员将海绵置于受试者头部与两侧线圈之间，以避免因受试者头部或身体移动产生运动伪影而降低图像质量。随后，叮嘱受试者闭眼、在扫描进行过程中尽量保持身体不动。

1.3 基于体素的形态学分析方法 本文对受试者的脑结构数据利用SPM（http://www.fil.ion.ucl.ac.uk/spm）软件及其VBM8 插件包进行全脑形态学分析。基于体素的形态学测量分析（voxel based morphometry,VBM）是一种在体素水平对全脑结构形态进行分析的数据驱动处理方法，其能定量评估脑局部灰质、白质和脑脊液密度或体积上的改变，从而精确测量局部脑组织形态学变化。该方法核心是通过图像分割技术提取全脑灰质、白质和脑脊液组织成分，再利用Dartel 方法将其配准至同一维度空间的标准模板，然后对全脑数据进行空间平滑，最后获得逐个体素的脑灰质体积（gray matter volume,GMV）指标（图1）。

图1 VBM 灰质体积计算示意图Fig 1 Schematic diagram of gray matter volume calculation(VBM)

1.4 基于KL 散度的个体化脑结构网络构建 本文应用GRETNA 软件（http://www.nitrc.org/projects/gretna）并结合上述脑灰质体积指标构建个体化全脑的结构网络（图2）。脑结构网络由节点（脑区）和边（脑区间形态特征的相似性）构成。

首先，网络节点的定义方式为：使用AAL（Anatomical Automatic Labeling） 模板将大脑分为90 个皮层和皮层下感兴趣区（region of interest,ROI），每个ROI 即代表一个网络节点。其次，网络边的定义为计算任意两个脑区间灰质体积概率密度分布相似度。其主要过程为：首先，按照上述网络节点定义提取各ROI 内所有体素GMV 值；随后，通过核函数估计方法获得各ROI 内所有体素GMV 值的概率密度分布（采样点个数为212=4 096）；最后，通过KL 散度衡量任意两个ROI 的GMV 概率密度分布之间相似度[9]，其定义如下。

如果随机变量X 的两种概率密度函数分别为P 和Q，则两者之间的KL 散度为：

图2 基于KL 散度的个体化脑形态相似度脑结构网络构建Fig 2 Brain structural network construction of individualized brain morphology similarity based on KL divergence

其中，KL 散度特性是KL（P，Q）≥0，当且仅当P=Q 时，KL（P，Q）=0。考虑这一特性，KL 散度可以在某种程度上反映两个概率密度函数间的相似性。但是，由于KL 散度不满足对称性要求，即KL（P，Q）≠KL（Q，P），所以任意两个脑区间的相似程度表示为：

通过指数变换将所有计算获得的任意两个脑区间的相似程度归一化，该值即表示两个脑区间KL 相似度（Kullback-Leibler similarity, KLS）系数，其公式表示为：

根据上述脑区之间相似度的定义，便可获得每位受试者90×90 的连接矩阵，即脑结构连接矩阵。根据不同的阈值选择，对每个脑结构矩阵进行阈值化（阈值的选择见下文），将其转化为二值矩阵，若两个脑区之间的结构连接大于设定的阈值，边就定义为1，两个脑区之间的结构连接小于设定的阈值，边就定义为0。

1.5 脑结构网络分析 本文采用多稀疏度阈值来对相似性矩阵进行阈值化，稀疏度是指一个网络内实际存在边连接数目与该网络内可能存在边连接总数间的比率。该方法保证所有受试者结构网络包含相同数目的边数，进而在相同网络结构条件下比较各网络拓扑指标的组间差异。因此，本文采用的稀疏度阈值范围是0.10～0.50，间隔为0.05。

在每个稀疏度阈值下，计算脑结构网络的全局指标，其主要包括5 个小世界属性指标[10]和2 个网络效率指标[11-12]。其中，小世界属性指标包括聚类系数Cp（反映网络的局部密度或专一化程度）、特征路径长度Lp（反映网络的平均连接或全局路径效率）、标准化聚类系数γ（真实网络和随机网络聚类系数的比值，γ>>1）、标准化特征路径长度λ（真实网络和随机网络特征路径长度的比值，λ≈1）、小世界属性σ=γ/λ（网络小世界属性的定量指标，σ>>1）[10]。另外，网络效率指标包括全局效率Eg（网络内并行信息传递的全局效率）和局部效率Eloc（网络的容错性）[11-12]。

进一步验证上述全局指标的可靠性，本文将上述真实网络的全局指标与通过一系列度分布相同的随机网络生成的全局指标进行若干次比较（例如，比较次数为1 000）判断真实网络的全局指标是否存在显著性以确保网络指标的可靠度。

1.6 统计学分析 为了评价两组受试者在年龄、性别和教育年限指标是否存在组间差异，本实验以组别作为固定因子，首先采用双样本t检验统计方法来确定教育年限和年龄两种因素是否具有组间差异；其次，通过χ2检验的统计方法评价两组性别因素是否存在组间差异。为了考查两组受试者脑结构网络的拓扑属性是否存在差异，本文用双样本t检验对每个稀疏度阈值下的全局网络指标（小世界属性和网络效率）进行比较。为了消除年龄、性别和教育年限因素对结果的影响，在所有分析中均把上述3个因素作为协变量进行回归。在显著性方面，本文使用的统计显著性为进行多重比较校正后阈值（P=0.05/9≈0.006）。另外，本文在统计分析过程中也选择适当宽松阈值进行组间比较，显著趋势阈值设为P＜0.05。

2 结果

2.1 人口统计学结果 本实验共纳入348 例试者，其中正常老龄组178 例，易转化AD 的MCI 高危患者组170 例，如表1 所示。统计结果表明，两组性别因素不存在显著性差异（χ2= 0.097，P=0.756）。两组间年龄和教育年限因素都不具有显著性差异（t=-1.126，P= 0.261；t=-1.021，P=0.308）。其中，正常老龄组的平均年龄为（75.43±5.00）岁，受教育程度为（16.08±2.88）年；易转化AD 的MCI 高危患者组的平均年龄为（74.68±7.25）岁，受教育程度为（15.77±2.85）年。

表1 人口统计学结果Tab 1 Demographic results

2.2 全局指标组间差异在定义的稀疏度阈值范围内，易转化AD 的MCI 患者组和正常老龄化对照组的脑结构网络与其对应的随机网络相比都具有较高的聚类系数（γ＞＞1）和近似相等的特征路径长度（λ≈1）。因此，结果表明两组受试者的脑结构网络都满足典型的小世界属性（σ＞＞1），如图3 所示。

虽然两组受试者的脑网络都满足小世界属性，但是通过双样本t检验分析发现两组受试者在小世界属性指标和网络效率上都存在显著性差异（图3、4）。与正常老龄对照组相比，在稀疏度阈值0.3～0.5 的范围内，易转化AD 的MCI 患者组的小世界性σ 显著减低，而在0.1～0.25 的范围内无明显差异（图3）。然而，患者组标准化的聚类系数γ 和特征路径长度λ 在全部稀疏度阈值的范围内均呈现显著减低（图3）。此外，易转化AD 的MCI 患者组的聚类系数Cp 和特征路径长度Lp 相比正常对照组呈现显著减低（图4）。与正常老龄对照组相比，在全部稀疏度阈值范围内，易转化AD 的MCI 患者组的全局效率Eg和局部效率Eloc 均呈现显著减低（图4）。

图3 不同稀疏度阈值下的小世界属性指标Fig 3 Small worldness index in the different sparsity thresholds

图4 不同稀疏度阈值下的全局指标Fig 4 Global index in the different sparsity thresholds

3 讨论

传统脑结构网络构建方法主要是利用脑结构（灰质）协共变网络分析方法进行构建。该方法是通过预先划定脑图谱模板，提取多名受试者结构像中每个脑区的形态测量特征（如GMV 值）；随后，通过上述模板依次计算每个脑区与全脑其它脑区GMV 值在跨被试间的共变关系来构建脑区间的联系[4-5，13]。本文采用的基于KL 散度的方法可以获得每名受试者的个体脑结构网络。前者主要依赖于大量受试者的脑结构数据，通过受试者组水平的分析方法进行脑网络构建，仅考虑受试者群体效应，而忽略受试者脑结构的个体差异。后者在个体水平获得结果能够集合多种统计方法（如方差分析，多元回归分析，双样本t检验，配对t检验等）并结合临床和认知行为解释大脑工作机制。这些为研究大脑结构的个体差异开辟新思路，特别是应用于识别每位患者的大脑异常情况。

本文运用基于KL 散度的方法构建个体脑结构网络，结合图论的方法以观察易转化为AD 的MCI患者的脑结构网络拓扑属性的改变。两个主要发现为：（1）基于KL 散度的个体化脑形态相似度网络的方法能够呈现人脑结构的拓扑属性，这为个体化精准诊疗提供新思路；（2）在全局水平，易转化AD 的MCI 患者组存在标准化聚类系数、标准化特征路径长度、小世界性、聚类系数、特征路径长度、全局效率和局部效率的显著减低，提示该疾病患者的脑结构网络向随机化转变。

虽然从生理学角度解释脑区之间的联系是相当复杂的，但是近期的研究试图推测脑区之间形态学相似性的理论基础。其中，对相似性的一种解释来源于大脑轴突理论[14]。该理论表明，在发育过程中由于脑区间连通性受到外力作用，使得脑区间产生相似的形态学特性。另一种解释表明，大脑的形成受到先天发育[15-16]和后天学习及记忆[17]等可塑性影响，使得各个脑区之间具有相似形态分布，进而反映了脑区之间的相互协调，以促使大脑完成各种认知和感知行为。综上所述，各个脑区之间形态学相似性关系可能被认为表征各个神经元或神经元聚类之间组成回路的结构特征。另外，各个脑区之间形态学相似性关系也可能说明相似程度主要取决于各个脑区内灰质细胞的种类类别（如，神经元，少突胶质细胞和星形胶质细胞等）和基因表达的分布情况[18]。然而，这些解释只是推测，尚未得到生物学研究的证明。因此，各个脑区之间形态学相似性关系的生物学基础将成为未来研究的重要课题。

人脑网络的主要特征是具有小世界网络属性，即具有较高的聚类系数和较短的路径长度，并兼具规则网络和随机网络的各自优势。与其相应，规则网络路径长度较长，具有较高的聚类系数；随机网络具有较短的路径长度，但是聚类系数较低。和正常老龄对照组相比，易转化AD 的MCI 患者组在部分稀疏度条件下小世界属性（σ）存在显著性减低，但是在全部稀疏度范围内标准化聚类系数（γ）和特征路径长度（λ）都明显小于正常对照组。另外，易转化AD 的MCI 患者组聚类系数、特征路径长度、局部效率和全局效率均存在显著性减低，提示其网络结构可能向随机化转变，这与早期关于AD 功能网络和解剖网络的研究表明AD 脑网络出现随机化转变的结果相一致[5,19]。全局效率的减少或者长距离连接的增加都会引起网络随机化。本文结果证明易转化AD 的MCI 患者的脑网络随机化是由上述因素引起的。因此，基于KL 散度的个体化脑形态相似度网络并结合图论有可能作为潜在的生物学标记用来监测易转化AD 的MCI 患者的疾病过程。