姚楠 苏春旺 李尤君 王珏 周昌松 黄子罡‡
1)(西安理工大学理学院,西安 710048)
2)(西安交通大学生命学院,生物医学信息工程教育部重点实验室,西安 710049)
3)(兰州大学物理科学与技术学院,兰州 730000)
4)(香港浸会大学物理系,非线性研究中心,计算及理论研究所,香港)
大脑具有自适应、自组织、多稳态等重要特征,是典型的复杂系统.人脑在静息态下的关键功能子网络——默认模式网络(DMN)的激活处于多状态间持续跳转的非平衡过程,揭示该过程背后的动力学机制具有重要的科学意义和临床应用前景.本文基于功能磁共振获得的血氧水平依赖(BOLD)信号,建立了DMN吸引子跳转非平衡过程的能量图景、吸引子非联通图、跳转关系网络等;以高级视觉皮层和听觉等皮层活动为例,通过对应激活DMN状态空间的分布,以及XGBoost、深度神经网络等算法验证了DMN状态变化与外部脑区状态的密切依赖关系;通过偏相关、收敛交叉映射等方法分析了DMN内各个脑区之间的相互作用.本文结果有助于理解静息态下大脑内在非平衡过程的动力学机制,以及从动力学的角度探索具有临床意义的脑功能障碍生物标志物.
近年来,非侵入式脑功能成像技术的发展促使脑科学研究进入了新的阶段[1].其中,fMRI(functional magnetic resonance imaging)可通过采集大脑的低频(< 0.1 Hz)血氧水平依赖(blood oxygen level dependent,BOLD)信号,来反映各个脑区在特定任务态或静息态下的动态激活状态.因无创和高空间分辨率等优势,fMRI被广泛应用于脑的工作机制和脑功能障碍的研究[2-5].基于fMRI对大脑在任务态和静息态等状态下激活模式的研究发现了大脑的多种功能子网络,包括默认模式网络(default mode network,DMN)、突显网络(salience network,SN)、执行控制网络(executive control network,ECN)、注意网络(attention network,AN)、视觉网络(visual network,VN)、听觉网络(auditory network,Aud)等[6].
DMN的正常运转对大脑功能的维持具有重要的意义.2001年,Raichle等[7]提出了大脑的“默认模式”(default mode)的概念,其对应的脑区在任务态或受到明显外部刺激时被抑制,而在静息态时恢复并维持激活,“默认模式”被认为是大脑神经系统的基线状态.2003年,DMN被明确提出[8],它是维持该基线状态最核心的功能子网络,DMN网络内不同脑区的静息态神经活动具有特定的时间相干性.此外,静息态下脑区之间的功能连接也被发现与对应的白质纤维束物理投射关系具有一定程度上的一致性[9].DMN的主要特征是在静息态下呈现出连续的高水平激活,而在任务态时受到抑制[10],已经证实它与维持觉醒和监测外部刺激等功能有关[11].在DMN的概念提出后,更多与大脑基线水平,即默认模式相关的脑区及脑区的具体功能被逐渐发现(见图1[12]).目前认为DMN主要包括: 腹内侧前额叶皮质(Ventral Medial Prefrontal Cortex,VMPC)、背内侧前额叶皮质(Dorsal Medial Prefrontal Cortex,DMPC)、后扣带回皮质(Posterior Cingulate Cortex,PCC)、楔前叶(Precuneus)和横向顶叶皮质(Lateral Parietal Cortex,LPC).
图1 大脑默认模式网络(DMN)[12]Fig.1.DMNof the brain[12].
DMN是大脑基线状态的核心功能网络,一些常见的脑功能障碍往往同时对应DMN激活模式及功能连接的异常,这也使得DMN活动的异常指标成为进行脑疾病的诊断和评估的潜在生物标志物[13-15],如与正常衰老相比,早期阿尔兹海默症(AD)患者的DMN中就存在显著的异常激活[16],轻度认知障碍(ADHD)患者在任务态下表现出DMN去激活率明显降低[17]和DMN内网络相关性降低[18],精神分裂症的阳性症状严重程度直接与DMN激活模式及功能连接的异常相关[19,20],脑外伤导致的轻度脑损伤患者的DMN则无法在任务态实现去激活,显著影响任务表现[21].
近年来,对静息态脑活动的研究发现,DMN单个脑区的激活看似随机,但从网络的层次却有规律可循.静息态下DMN的状态存在更细致的子状态,对应于不同的激活模式,系统在多种激活模式间持续跳转,是典型的非平衡过程,其表现出的规律及背后的机理适合运用统计物理和非线性动力学的方法进行研究.然而,在这一方向上更进一步的深入研究到目前为止还尚未见到进展,例如,静息态DMN的子状态是否和其他脑区的功能活动有对应关系等,需要做进一步的考察.此外,由于脑功能障碍对应的DMN激活模式及跳转行为与健康状态存在显著差异[13,19,20,22],对静息态下DMN所处非平衡过程的研究具有重要的科学意义和临床应用前景,能够帮助我们理解人脑静息态的动力学特征和工作机制,从功能连接的角度揭示脑功能障碍的发病机理,并为相关疾病的预防、诊断和治疗提供启发[22,23,24].
鉴于上述原因,本文将基于静息态DMN的BOLD数据,围绕DMN的动力学行为及其背后的机理展开研究,所采用的方法主要包括能量图景构造[25,26]、共变行为分析、基于压缩感知的偏相关分析[27]与收敛交叉映射方法[28]等.研究的问题包括:DMN激活行为的非平衡过程的主要特征,DMN与外部其他脑区活动的交互,DMN内部各脑区之间的相互作用,以及对其非平衡过程动力学机理的探讨.
2.1.1 数据及预处理
本文的研究采用静息态功能磁共振成像(rsfMRI)数据来源于人类连接组计划(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.
表1 默认模式网络ROI的名称、标签、Brodmann分区(BA)编号及位置坐标信息Table 1.Name,label,Brodmann area(BA)number,and location information of ROIs belonging to the default mode network.
2.1.2 DMN的能量图景构造方法
以背侧默认模式网络(dDMN,为DMN的核心区域)作为代表,将其9个ROI的BOLD信号根据某阈值(一般取均值线或以均值线为中心的窄带)二值化为0或1序列,分别对应于不活跃状态和活跃状态,如图2所示,编号为i的ROI的状态为 σi=0,1(i=1,2,···,9),那么 t时刻 dDMN状态可表示为 σ =[σ1,σ2,···,σ9] .系统状态s所对应的“能量”定义为
其中J ∈R9×9和h ∈R9是模型的两组参数,Jii=0(i=1,2,···,9).以上能量函数给出系统状态的玻尔兹曼分布:
配分函数为
在真实rs-fMRI数据dDMN脑区激活占比的约束条件下,通过梯度上升算法,不断调节(即训练)模型的两个参数J和h,使得玻尔兹曼分布的熵
局部极小态及能量壁垒的计算: 将dDMN状态矢量中只有一位存在差异的两个状态视为相邻状态,因此每个dDMN状态都有9个相邻状态.若某状态的能量小于所有9个相邻状态,则称为局部极小态.
局部极小态的非联通图构建: 1)首先建立状态矢量构成的超晶格,其中每个状态都与其相邻状态相连;2)设定能量阈值,Eth,它等于当前所有状态中最大的能量;3)移除所有能量 ≥ Eth的状态及其连边;4)检测超晶格中是否每个局部极小态都有至少一条路径相互联通;5)重复步骤3)和4),将Eth设为剩余格点中的能量最大值,直到所有的局部极小态都互不联通.记录下每两个局部极小态第一次不连通时所设定的Eth,作为两者之间的势垒高度,即可构建局部极小态的非联通图.
图2 DMN能量图景构造方法示意图Fig.2.Schematic diagram of the construction of DMN energy landscape.
2.1.3 基于压缩感知的偏相关方法
基于脑区BOLD信号的相关性分析是构建脑功能网络的常用方法.但相关性分析所得两变量(如DMN内部某两个ROI)间关系无法排除来自其他变量(即除这两个ROI之外的DMN内部其他ROI)的干扰,而偏相关分析则可以在一定程度上去除其他变量的影响.压缩感知(包括LASSO)利用系统稀疏性的先验知识辅助求解线性回归问题,在小样本数据集的情况下优势更加突出,甚至可解线性欠定方程组.脑功能网络具有稀疏性这一基本特征,而且核磁共振数据的时间序列长度一般都比较短,利用压缩感知辅助进行偏相关分析构建脑功能网络,不仅可以凸显网络中的重要边、自动剔除弱的伪连边,而且还可以对基于短时间序列的相关性(或线性回归)计算的准确性起到一定的补偿.以下将简要介绍基于压缩感知的偏相关方法[27].
假设 f1,f2,···,fp是所测p个对象的数据矢量(时间序列),长度为 n,且已经被中心化(centered,即其中 1n∈ Rn×1为所有元素均为1的n维列矢量)、标准化(normalized,即那么协方差矩阵为
协方差矩阵的逆矩阵(称为精度矩阵,precision matrix)为
则fi与fj之间的相关系数和偏相关系数分别为
如果对fi作如下线性回归:
那么回归系数αij就等于fi与fj之间的相关系数ρij.而如果对fi做如下线性回归:
那么回归系数βij和fi与fj之间偏相关系数θij的关系为
其中var(ϵ)表示回归残差ϵ的方差.
在回归系数矢量β具有比较强的先验稀疏性(在脑功能网络中这是一个合理的假设)的情况下,可被如下LASSO回归的方法很好地估计:
在通过上述方法得到β之后,即可计算偏相关系数,
其中
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的结构得到预测序列,求得与实际序列Y的皮尔森相关系数ρY|MX,此系数即是变量Y对变量X影响程度的度量;同样地,ρX|MY是变量X对变量Y影响程度的度量.X与Y之间的因果方向根据R=ρX|MY- ρY|MX的符号进行判断(R > 0则表明X对Y的动力学影响更强,R < 0则反之).
2.2.1 DMN静息态非平衡过程的动力学行为分析
2.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(d)—图3(i)展示了55个被试能量图景的统计分布特征.可以看到: 吸引子数目存在明显的峰值(图3(d),峰值位于4附近),表明静息态下吸引子典型数目的普遍存在.吸引域占据率(即被访问次数)呈双峰分布(图3(e)),表明静息态下DMN的状态多处在其中某些吸引域而很少在另外的吸引域,存在两极分化现象.类似地,势阱深度(即势阱最低点到最近鞍点的能量差)也呈双峰分布(图3(f)),势阱深度和吸引域占据率的正相关关系如图3(g)所示.图3(h)给出了吸引域被访问次数O和吸引域面积A(即吸引域内包含的状态总数)之间具有很好的线性关系 O=kA(其中k为常数),而具有较大面积的吸引域也往往具有较大的势阱深度,其具体关系曲线见图3(i),这与图3(g)中吸引域被访问次数和势阱深度之间的正相关对应.结合DMN能量图景结构及跳转过程可以得到以下启发.
图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)能量壁垒和吸引域面积的关系Fig.3.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.
a)不同被试的能量图景吸引子数目以及吸引域间跳转行为存在差异,分别反映了其在静息态下脑活动的丰富程度和活跃程度的差异.这可能与不同被试的认知及思维活动特征差异有关,也可能与被试在实验测量前的短期经历有关.已有研究发现静息态数据采集之前被试所处环境对其静息水平存在遗留影响.可以推测,被试在进入静息态之前若经历了较为丰富的刺激,所得到的DMN静息态能量图景应具有更多的吸引子状态和更频繁的跳转.该推测还有待进一步的实验验证.
此外,能量图景的吸引子势能、位置关系和吸引域大小等信息可提示潜在的激活模式间跳转方式.如被试1的数据中存在跨越高势垒(见图3(a))的状态2—状态5间频繁交替的模式(状态5具有较大的吸引域),跨越低势垒的状态3—状态4间频繁交替的相对局部的模式(其状态较为接近),以及由上述2—5相到3—4相过渡的状态3—状态5交替的中间模式.据此可以了解某被试在静息态下其DMN状态的特有跳转模式,由于DMN的各个脑区激活行为与外部脑区的神经活动状态存在关联,可以根据该跳转模式理解被试在静息态下的习惯思维活动.
b)目前基于 BOLD数据对 fingerprint(指纹)指标,即稳定的个体间差异指标的探讨是一个重要的前沿科学问题[29].从上述分析可以看到,不同被试的吸引子个数、状态空间轨迹构型、跳转方式等均具有特异性,其中对特定个体相对稳定同时在不同个体之间差异较大的动力学属性可以作为潜在的动力学指纹,并运用于对被试的个体属性及脑功能障碍的识别.
2.2.1.2 DMN与外部其他脑区活动的交互
以上介绍了静息态下DMN状态在不同吸引子间跳转的行为,该非平衡过程的机理解释需要从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三种方法及两种训练方式Fig.4.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.
图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轨道显著的个体间差异所影响.
图5 DMN状态与视听觉皮层状态的对应现象(a),(b)根据吸引子的吸引域标签着色的DMN状态轨道;(c),(d)根据(hV,Aud)标签序列着色的DMN状态轨道;(a),(c)被试1的结果;(b)(d)被试2的结果Fig.5.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特定形式的抑制作用.
2.2.2 静息态DMN内部功能网络分析
为了分析不同被试在静息态DMN内部功能网络(这里包括关联网络和因果性网络)方面的共性特征,我们计算了平均的结果,得到了一些在静息态下DMN内部普遍存在的功能性结构特征.
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(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为反相关.
2.2.2.2 静息态DMN内部的因果网络
图6(c)为CCM方法中得到的i脑区对j脑区的影响程度ρi|Mj,该矩阵为近似对称的非对称矩阵.可以看到,其中较强的影响程度对应的脑区在偏相关分析中(图6(a))也被发现具有较强的关联性.图6(d)给出了CCM方法得到的动力学因果性网络,可以明显地看到DMN的子网络结构,其中vDMN,pDMN相对分离,dDMN内部D1-D6相对一致,D7(丘脑)、D8和D9(左右海马旁回)则相对独立.CCM方法得到的结果R可对应于DMN内部区域间的信息流,见图6(e).在静息态下,丘脑(D7)和左右海马旁回(D8,D9)处于DMN内部信息流的上游,dDMN中的D1—D6和pDMN是下游子网络;vDMN从丘脑和左右海马旁回接收信息,并向pDMN和dDMN的主区域(D1—D6)发送信息.在dDMN中,D1(内侧额叶)相对特殊,它对pDMN和vDMN中的多个脑区都有较强的动力学影响;在vDMN中,V3和V8(左右嘴峰,Culmen)也具有广泛的对其他脑区的影响.
研究表明,大脑DMN在静息态下出现多种激活模式切换的行为,但从复杂系统非平衡动力学的角度研究静息态DMN激活行为的工作还比较少见.本文基于静息态fMRI数据分析了静息态DMN的动力学行为.首先通过构造静息态DMN能量图景,观察了不同被试在静息态的DMN状态非平衡演化过程.为了理解该非平衡过程受到外部影响的机制,我们以高级视觉皮层和听觉皮层为例,通过DMN状态空间中的分布特征以及XGBoost、深度神经网络等机器学习算法,明确了外部脑区活动与DMN状态跳转的密切关系.这为深入理解DMN的动力学行为启发了新的思路.文中图3与图4的结果也验证了对DMN的传统认识,即DMN作为维持觉醒功能的网络会因外部刺激而发生状态调节.并且,我们的结果从能量图景与外部脑区激活模式标签对照的新视角提供了静息态下DMN的动力学状态与外部脑区(视、听觉皮层)状态之间关联的证据.
另外,本文采用偏相关、收敛交叉映射等方法计算了DMN内部23个ROI的功能网络,包括关联性网络和因果性网络,识别出了部分主要功能性连接,以及DMN内部的子网络结构和动力学影响的核心区域(见图6).虽然DMN的动力学行为具有较大的个体间差异,但组平均的关联网络能够反映出比较稳定的强相关边;因果性网络给出了DMN的子网络结构,该结构可以得到关联性网络的佐证.这些结果印证了前人对于DMN主要内部结构的研究结论(例如DMN分为dDMN,vDMN和pDMN三大区域,其各自具有密集的内部连接而彼此之间连接相对较少,丘脑和海马旁回则是相对独立的功能区等).同时,也给出了一些新的启发,例如根据因果性网络平均结果给出的DMN内部ROI间信息流,我们可以得到这样的结论: 在静息态下,1)丘脑和海马旁回是DMN其他区域的驱动者,即信息是从丘脑和海马旁回流向DMN其他区域;2)嘴峰(Culmen,相应的ROIs为vDMN3和vDMN8)是DMN内很多其他区域的驱动者(甚至处于丘脑和海马旁回的上游).
上述围绕DMN内部功能连接以及与外部脑区共变特征的探讨,增加了我们对DMN动力学行为的理解.但由于DMN的吸引子状态激活与视听觉、突显网络等功能网络活动的相互作用具体细节机制尚不完全清楚,目前还无法把外部脑区对DMN活动非平衡过程的影响形式化为特定的方程,这还有待神经影像学与复杂系统等领域的联合探索.本文基于静息态DMN状态能量图景的动力学分析方法具有潜在的应用价值,例如: 1)已有研究发现静息态数据采集之前,被试所处环境对其静息水平存在遗留影响(例如,被试刚刚经历了嘈杂环境或完成任务态测试,与充分休息后所观察到的静息态存在一定差异),那么,特定被试在不同“静息水平”下的DMN状态分布及切换行为就可以根据状态能量图景中各个吸引子的能量大小以及吸引域大小进行分析.例如,“弱静息”条件下(高“温度”),DMN可能更多地在具有较大吸引域的亚稳态附近漫游,而“强静息”条件下(低“温度”)则陷入具有较低“能量”的亚稳态附近,即使其吸引域较小;2)根据静息态DMN状态能量图景中各个吸引子相对位置以及吸引域构型特征,结合DMN状态与外部脑区激活的对应关系,预测被试在静息态下的特征“心理活动轨迹”.
感谢北京师范大学王文旭教授,兰州大学俞连春、黄亮教授,陕西师范大学王圣军教授的有益讨论,以及西安交通大学黄文敏、曹玲灿、李蒙在数据处理方面的协助.