全文HTML
--> --> --> -->2.1.方 法
32.1.1.数据及预处理
本文的研究采用静息态功能磁共振成像(rs-fMRI)数据来源于人类连接组计划(human connectome project, HCP)网站的公开数据集, fMRI空间分辨率为2 mm × 2 mm × 2 mm (单个体素尺寸), 时间分辨率为0.72 s (采样时间间隔), 包括55个被试, 每个被试的采样点数为T = 1200, 采样持续时间为0.72 s × 1200 = 14 min. 原始影像数据的预处理流程为: 功能与结构图像配准, 配准到标准空间, 去除白质和脑脊液信号, 头动校正, 带通滤波(0.01—0.10 Hz), 空间高斯平滑, 去漂移等. 经上述预处理后, 再根据脑区划分模板将默认模式网络相关的区域划分为23个感兴趣区 (regions of interest, ROIs), 并通过对ROI内所有体素时间序列求平均, 得到每个ROI的BOLD(血氧水平依赖)信号时间序列. 23个ROI的相关基本信息见表1.ROI | Label | L/R | BA | X | Y | Z |
Posterior Cingulate | vDMN_1 | L | 31 | –12 | –62 | 10 |
Middle Frontal Gyrus | vDMN_2 | L | 10 | –27 | –6 | 59 |
Culmen | vDMN_3 | L | 37 | –30 | –39 | –20 |
Superior Occipital Gyrus | vDMN_4 | L | 19 | –36 | –88 | 28 |
Posterior Cingulate Gyrus | vDMN_5 | R | 31 | 15 | –56 | 13 |
Precuneus | vDMN_6 | 7 | –6 | –61 | 56 | |
Middle Frontal Gyrus | vDMN_7 | R | 10 | 24 | 26 | 47 |
Culmen | vDMN_8 | R | 37 | 27 | –33 | –23 |
Angular Gyrus | vDMN_9 | R | 39 | 43 | –79 | 28 |
Cerebellum | vDMN_10 | R | 12 | –47 | –63 | |
Ventral Posterior Cingulate Gyrus | pDMN_1 | 23 | 0 | –35 | 28 | |
Precuneus | pDMN_2 | 7 | 0 | –76 | 38 | |
Inferior Parietal Lobule | pDMN_3 | L | 40 | –39 | –64 | 46 |
Inferior Parietal Lobule | pDMN_4 | R | 40 | 39 | –64 | 46 |
Middle Frontal Gyrus | dDMN_1 | 9 | 0 | 49 | 12 | |
Angular Gyrus | dDMN_2 | L | 39 | –48 | –73 | 32 |
Superior Frontal Gyrus | dDMN_3 | R | 6 | 18 | 38 | 51 |
Dorsal Posterior Cingulate Gyrus | dDMN_4 | 31 | 0 | –57 | 30 | |
Ventral Anterior Cingulate Gyrus | dDMN_5 | 24 | 0 | –17 | 35 | |
Angular Gyrus | dDMN_6 | R | 39 | 48 | –66 | 29 |
Thalamus | dDMN_7 | –6 | –6 | 3 | ||
Parahippocampal Gyrus | dDMN_8 | L | 36 | –24 | –37 | –9 |
Parahippocampal Gyrus | dDMN_9 | R | 36 | 24 | –21 | –23 |
表1默认模式网络ROI的名称、标签、Brodmann分区(BA)编号及位置坐标信息
Table1.Name, label, Brodmann area (BA) number, and location information of ROIs belonging to the default mode network
3
2.1.2.DMN的能量图景构造方法
以背侧默认模式网络(dDMN, 为DMN的核心区域)作为代表, 将其9个ROI的BOLD信号根据某阈值(一般取均值线或以均值线为中心的窄带)二值化为0或1序列, 分别对应于不活跃状态和活跃状态, 如图2所示, 编号为i的ROI的状态为σi = 0, 1(图 2 DMN能量图景构造方法示意图
Figure2. Schematic diagram of the construction of DMN energy landscape.
局部极小态及能量壁垒的计算: 将dDMN状态矢量中只有一位存在差异的两个状态视为相邻状态, 因此每个dDMN状态都有9个相邻状态. 若某状态的能量小于所有9个相邻状态, 则称为局部极小态.
局部极小态的非联通图构建: 1)首先建立状态矢量构成的超晶格, 其中每个状态都与其相邻状态相连; 2)设定能量阈值, Eth, 它等于当前所有状态中最大的能量; 3)移除所有能量 ≥ Eth的状态及其连边; 4)检测超晶格中是否每个局部极小态都有至少一条路径相互联通; 5)重复步骤3)和4), 将Eth设为剩余格点中的能量最大值, 直到所有的局部极小态都互不联通. 记录下每两个局部极小态第一次不连通时所设定的Eth, 作为两者之间的势垒高度, 即可构建局部极小态的非联通图.
3
2.1.3.基于压缩感知的偏相关方法
基于脑区BOLD信号的相关性分析是构建脑功能网络的常用方法. 但相关性分析所得两变量(如DMN内部某两个ROI)间关系无法排除来自其他变量(即除这两个ROI之外的DMN内部其他ROI)的干扰, 而偏相关分析则可以在一定程度上去除其他变量的影响. 压缩感知(包括LASSO)利用系统稀疏性的先验知识辅助求解线性回归问题, 在小样本数据集的情况下优势更加突出, 甚至可解线性欠定方程组. 脑功能网络具有稀疏性这一基本特征, 而且核磁共振数据的时间序列长度一般都比较短, 利用压缩感知辅助进行偏相关分析构建脑功能网络, 不仅可以凸显网络中的重要边、自动剔除弱的伪连边, 而且还可以对基于短时间序列的相关性(或线性回归)计算的准确性起到一定的补偿. 以下将简要介绍基于压缩感知的偏相关方法[27].假设
在回归系数矢量β具有比较强的先验稀疏性(在脑功能网络中这是一个合理的假设)的情况下, 可被如下LASSO回归的方法很好地估计:
3
2.1.4.收敛交叉映射方法(CCM)
鉴于大脑活动的典型非线性特征及fMRI所得BOLD数据采集时长有限两方面因素, 我们还将采用收敛交叉映射法(convergent cross mapping, CCM)[28]进行动力学因果性分析. CCM基于泰肯嵌入定理(Taken's embedding theorem), 通过相空间重构的方式对变量间的动力学因果性进行推断. 相较于目前广泛应用的格兰杰因果性分析(Granger causality, GC)和传递熵(transfer entropy, TE)等方法, CCM的优势在于更加适合非线性系统的处理, 并且对数据量的要求相对较低.通过CCM方法对两个脑区BOLD时间序列X和Y进行因果性分析, 可以得到脑区间有向的动力学影响强度(directed dynamical influence), 进而构建DMN与其他脑区相互作用的脑网络. 基本步骤包括: 由X和Y时间序列分别进行相空间重构得到流形MX与MY; 根据MX的结构得到预测序列
2
2.2.结 果
32.2.1.DMN静息态非平衡过程的动力学行为分析
32.2.1.1.DMN的吸引子间跳转行为特征
基于前面所述方法, 可由每个被试的数据计算得到dDMN能量图景. 这里, 高维超晶格上的能量图景的主要特征可通过非连通图(disconnectivity graph)显示, 包括吸引子能量大小及相邻最低势垒连通关系等信息, 见图3(a). 非连通图的构造是在训练得到系统能量函数表达式后, 通过搜索能量极小值, 分析各极小值对应的系统状态吸引子的位置关系, 最终保留所有吸引子及它们之间的最低翻越势垒构成的树状结构图. 图3(a)—图3(c)上下两行分别对应于被试1和被试2的结果. 以能量图景吸引域为标签对DMN活动时间序列进行着色可对其在状态空间中的分布进行观察, 见图3(b), 其中左右图为不同时刻的状态激活情况(为观察方便计, 这里使用了对同一状态轨道不同视角的展示, 即左右图观察角度不同). 三维投影空间根据23维数据主成分分析(PCA)得到的前三个主成分方向进行构造. 图3(c)给出了吸引子间状态跳转关系网络, 其中颜色与图3(a)中的吸引子一致, 节点大小等于DMN在该吸引域内出现的次数, 节点间有向边的权重表示对应方向的跳转频次, 总时长T = 1200. 比较图3(a)—图3(c)上下两行对应的两个不同被试的结果可以看到, 不同被试在静息态下其默认模式网络激活行为的动力学特征存在明显差异, 例如被试1具有较多的吸引子状态, 且存在更频繁的吸引子间切换行为, 而被试2的吸引子个数较少, 且跳转行为主要表现为局部切换(例如state 1和state 2之间, state 3和state 4之间), 这表明静息态下被试1的思维活动更加丰富、活跃, 而被试2的思维活动则相对单一且不易变化.图 3 默认模式网络(DMN)激活行为的动力学特征 (a) dDMN状态能量图景中吸引子状态的非联通图(disconnectivity graph); (b)根据吸引子吸引域标签进行着色的DMN状态轨迹在3维空间的投影, 左右两图显示了同一被试轨迹不同角度的观察; (c)吸引子间跳转网络, 节点大小表示对应吸引域的访问次数, 点间有向边的权重表示对应的单向跳转频次, 总时长T = 1200. (a)?(c)中上下两行分别对应于被试1和被试2的结果; (d)?(i)为55个被试能量图景的统计分布特征: (d)吸引子数目的统计分布; (e)吸引域占据率(访问次数)的统计分布; (f)能量壁垒(势阱深度)的统计分布; (g)吸引域占据率与能量壁垒的关系; (h) 吸引域占据率与吸引域面积的关系; (i) 能量壁垒和吸引域面积的关系
Figure3. Dynamic characteristics of the resting state default mode network (DMN): (a) Disconnectivity graph of the energy landscape for the dDMN; (b) projections of the DMN state trajectories, colored according to the labels of attraction basins, in a 3-dimensional space. The left and right panels show the trajectory of the same subject from different viewing angles, respectively; (c) attractors transition networks. The nodal size indicates the visiting frequency in the corresponding attracting domain, and the weight of directed edge indicates the corresponding directed switching frequency, with the recording time T = 1200. The upper and lower rows in the panels (a)?(c) correspond to the results of test subjects 1 and 2, respectively. Panels (d)?(i) are the statistical results of energy landscapes for 55 subjects: (d) The distribution of attractor numbers; (e) the distribution of basin occupations; (f) the distribution of energy barriers; (g) dependence of basin occupation and energy barrier; (h) dependence of basin occupation and area; (i) dependence of energy barrier and area.
图3(d)—图3(i)展示了55个被试能量图景的统计分布特征. 可以看到: 吸引子数目存在明显的峰值(图3(d), 峰值位于4附近), 表明静息态下吸引子典型数目的普遍存在. 吸引域占据率(即被访问次数)呈双峰分布(图3(e)), 表明静息态下DMN的状态多处在其中某些吸引域而很少在另外的吸引域, 存在两极分化现象. 类似地, 势阱深度(即势阱最低点到最近鞍点的能量差)也呈双峰分布(图3(f)), 势阱深度和吸引域占据率的正相关关系如图3(g)所示. 图3(h)给出了吸引域被访问次数O和吸引域面积A(即吸引域内包含的状态总数)之间具有很好的线性关系
a)不同被试的能量图景吸引子数目以及吸引域间跳转行为存在差异, 分别反映了其在静息态下脑活动的丰富程度和活跃程度的差异. 这可能与不同被试的认知及思维活动特征差异有关, 也可能与被试在实验测量前的短期经历有关. 已有研究发现静息态数据采集之前被试所处环境对其静息水平存在遗留影响. 可以推测, 被试在进入静息态之前若经历了较为丰富的刺激, 所得到的DMN静息态能量图景应具有更多的吸引子状态和更频繁的跳转. 该推测还有待进一步的实验验证.
此外, 能量图景的吸引子势能、位置关系和吸引域大小等信息可提示潜在的激活模式间跳转方式. 如被试1的数据中存在跨越高势垒(见图3(a))的状态2—状态5间频繁交替的模式(状态5具有较大的吸引域), 跨越低势垒的状态3—状态4间频繁交替的相对局部的模式(其状态较为接近), 以及由上述2—5相到3—4相过渡的状态3—状态5交替的中间模式. 据此可以了解某被试在静息态下其DMN状态的特有跳转模式, 由于DMN的各个脑区激活行为与外部脑区的神经活动状态存在关联, 可以根据该跳转模式理解被试在静息态下的习惯思维活动.
b)目前基于BOLD数据对fingerprint(指纹)指标, 即稳定的个体间差异指标的探讨是一个重要的前沿科学问题[29]. 从上述分析可以看到, 不同被试的吸引子个数、状态空间轨迹构型、跳转方式等均具有特异性, 其中对特定个体相对稳定同时在不同个体之间差异较大的动力学属性可以作为潜在的动力学指纹, 并运用于对被试的个体属性及脑功能障碍的识别.
3
2.2.1.2.DMN与外部其他脑区活动的交互
以上介绍了静息态下DMN状态在不同吸引子间跳转的行为, 该非平衡过程的机理解释需要从DMN的内在动力学以及外部其他脑区的影响两方面考虑. DMN作为众多脑功能网络之一, 其状态会受到外部其他功能网络的影响. 这里我们首先以高级视觉皮层(hV)和听觉皮层(Aud)为例, 分析DMN与其他脑区激活行为的交互关系.图4(a)给出了高级视觉皮层(hV)和听觉皮层(Aud)的BOLD时间序列, 若进行二值化, 用1和0对应激活与不激活, 则时间序列可用4种(hV, Aud)激活状态表示, 即(1, 0), (0, 1), (0, 0)和(1, 1), 图中用4种颜色示意. 为了分析DMN与hV, Aud的共变行为, 可根据同时刻的(hV, Aud)标签对DMN状态轨道进行着色, 如图4(b), 并与图3(b)中根据吸引域标签染色的轨迹进行对比. 图5(a)和图5(b)分别给出了被试1和被试2根据吸引域标签进行着色的DMN状态轨迹在2维空间(由PCA前两个主成分对应的本征矢张成)的投影, 图5(c)和图5(d)给出了被试1和被试2根据(hV, Aud)标签着色的同样的轨迹投影. 从图5可见不同的(hV, Aud)标签在DMN的状态空间中对应特定的DMN子状态区域(即吸引子及其吸引域), 说明DMN的状态与(hV, Aud)的状态有较好的对应关系. 我们通过15个被试数据的组分析对该现象进行验证, 也可见一致的现象, 见图4(c), 其中二维空间是由(通过对15个被试数据进行组PCA得到的)前两个主成分PCA1和PCA2张成, 每个点表示某被试某一时刻的DMN状态投影. 所有被试的DMN均表现出了一致地与自身(hV, Aud)存在对应关系的激活行为, 这一规律没有被DMN轨道显著的个体间差异所影响.
图 4 默认模式网络(DMN)状态与其他脑区激活行为的共变 (a)对高级视觉皮层(hV)和听觉皮层(Aud)活动二值化, 得到4种激活状态的(hV, Aud)标签序列; (b)根据(hV, Aud)标签序列着色的DMN状态轨道; (c) DMN的23个脑区激活行为在PCA1-PAC2空间的投影; (d)对DMN的激活行为数据按照同时刻(hV, Aud)标签进行4分类训练所得ROC曲线, 其中包括XGBoost, DNN, Random Forest三种方法及两种训练方式
Figure4. Covariation of the default mode network (DMN) states with activity mode in other brain regions: (a) Binarizing the activities of high-level visual cortex (hV) and auditory cortex (Aud) to obtain (hV, Aud) label sequence; (b) DMN state orbit colored according to (hV, Aud) label sequence; (c) projection of the activity states of the 23 brain regions in the PCA1-PAC2 space; (d) the ROC curve obtained by performing 4-classes training on the (hV, Aud) labels of DMN activity data by XGBoost, DNN, Random Forest with two training schemes.
图 5 DMN状态与视听觉皮层状态的对应现象 (a), (b)根据吸引子的吸引域标签着色的DMN状态轨道; (c), (d) 根据(hV, Aud)标签序列着色的DMN状态轨道; (a), (c)被试1的结果; (b)(d)被试2的结果
Figure5. Correspondence between the states of DMN and (visual, auditory) cortexes: (a), (b) DMN state orbits colored according to the labels of attraction basin; (c), (d) DMN state orbits colored according to the labels of (hV, Aud) states. (a), (c) the results of subject 1; (b)(d) the results of subject 2.
以上是从状态空间的分布特征观察DMN与(hV, Aud)的状态对应关系, 我们还可通过XGBoost, DNN(deep neural network), Random Forest等机器学习算法对上述关系进行验证. 以(hV, Aud)的4个状态作为标签, 用15个被试DMN的23个脑区激活行为的1200个采样数据对模型进行4分类训练. 图4(d)给出了对应的ROC曲线, 横坐标为假阳性率(FPR), 纵坐标为真阳性率(TPR), 图中虚线表示通过跨被试的交叉验证, 即随机选择13个被试构成训练集, 剩余的2个被试作为测试集的结果. 图中实线表示通过十折交叉验证得到的结果, 具体处理方法是对某被试的1200个数据点随机拆分为训练集和测试集, 用十折交叉验证得到结果, 再对15个被试的FPR和TPR进行平均得到ROC曲线. 从ROC曲线可以看出: 1)根据某被试DMN的状态能够得到对其自身(hV, Aud)状态的较为准确的推断, 即DMN与hV, Aud的状态存在密切的依赖关系; 2)跨被试的验证结果劣于单个被试十折交叉验证结果, 是由于DMN轨迹的个体间差异等因素所致.
从上述DMN与(hV, Aud)的对应关系可知, 静息态下DMN在多吸引子间跳转的动力学过程与DMN外部脑区活动密切相关. 需要说明的是, 静息态并非对人脑状态的严格定义, 而是BOLD数据采集的实验条件, 在静息态下, 某些事件(也可能是自发涨落)导致视觉或听觉皮层的激活, 会伴随着SN等网络参与的对DMN特定形式的抑制作用.
3
2.2.2.静息态DMN内部功能网络分析
为了分析不同被试在静息态DMN内部功能网络(这里包括关联网络和因果性网络)方面的共性特征, 我们计算了平均的结果, 得到了一些在静息态下DMN内部普遍存在的功能性结构特征.3
2.2.2.1.静息态DMN内部的关联网络
基于偏相关计算得到DMN的23个ROI的关系矩阵及网络(见图6(a)和图6(b)). 可以看到, 网络中强的正相关性边包括vDMN (V1-V10), pDMN (P1-P4), dDMN的主区域(D1-D6)以及左右海马旁回(parahippocampal gyrus, 图中的D8和D9)各自的内部连边. 此外, 强的正相关性连边还包括D2 (左角回, angular gyrus)-P3 (左下顶叶, inferior parietal lobule)、D4 (左背侧后扣带回, dorsal posterior cingulate gyrus)-P2 (右楔前叶, precuneus)、D4-V1 (左后扣带回, posterior cingulate)、D4-V5 (右后扣带回, posterior cingulate)、P2-V6 (左楔前叶, precuneus)、D5 (右腹侧扣带回前部, ventral anterior cingulate gyrus)-P1 (左腹后扣带回, ventral posterior cingulate gyrus)、P4 (右下顶叶, inferior parietal lobule)-V9 (右角回, angular gyrus)等, 这些强的正相关点对多为空间近邻.图 6 DMN内部23个ROI间关系分析 (a), (b)基于压缩感知的偏相关矩阵及网络结构, 网络中紫色线表示正相关, 青色线表示负相关, 线的粗细表示相关性强弱; (c) CCM方法中得到的ROI之间的单向影响程度
Figure6. Relationships of the 23 ROIs within DMN: (a) Partial correlation matrix calculated based on compressive sensing; (b) the corresponding network structure, which includes the positive (magenta links) and negative (cyan links) correlations. The width of the line indicates the correlation strength; (c) the directed dynamical influence from ROI i to j,
图6(a)和图6(b)中也存在强的反相关性边, 如D2-P1, P1-V9, D1-P4, D6-V4. 对于反相关性边, 一种可能的原因是信息沿着正相关性边所连成的通路传递的过程中积累的延迟相位差导致的, 并且两节点的反相关程度与该两节点之间由正相关边所组成的最短路径长度有关[30], 随着最短路径长度的增加, 反相关程度增强, 继续增加又会使其减弱(相位引起的效应). 这里我们从DMN的ROI中也观察到了类似的情况, 如D6-V9, V9-V4均为正相关, V6-D4为反相关; P1-P4, P4-V9均为正相关, P1-V9为反相关.
3