静息态下大脑的动态模块化指纹
摘要：人脑是一个动态的模块化网络，可以分解为一系列模块，其活动随时间不断变化。静息状态下，在亚秒级的时间尺度上会出现几个脑网络，即静息态网络(RSNs)，并进行交互通信。本文尝试探究自发脑模块化的快速重塑及其与RSNs的关系。三个独立的健康受试者静息态数据集（N=568），对其使用脑电/脑磁图(EEG/MEG)来探究模块化脑网络的动态活动。本文证实了RSNs的存在，且其中一些网络存在分裂现象，尤其是默认模式网络、视觉、颞区和背侧注意力网络。本文也证明了心理意象中的个体间差异与特定模块的时间特征有关，尤其是视觉网络。综上所述，本文的研究结果表明大规模电生理网络在静息态时具有依赖模块化的动态指纹。
1. 引言
自发脑活动在从亚秒到年的多个时间尺度上不断变化。这些波动涉及一系列网络，称为静息态网络(RSNs)。大多研究认为自发脑活动可以用一系列时空网络模式来解释。主要结合降维算法（如K均值聚类）、主成分分析、正交连通性因子分解或基于模型的方法（如隐马尔可夫模型）。从这些快速动态分析中得出的特征也可作为某些脑疾病和行为特征的潜在神经标志物。
有研究表明人脑是一个模块化网络，可划分为“模块”（也称为社区或聚类），即大脑内部连接紧密但与其他区域连接不紧密的区域。虽然有几项研究探究了大脑模块化网络的时间依赖性，但在与快速认知相关的时间尺度上，对自发模块化脑网络的快速重塑及其与RSNs的关系的研究是很少的。为了精确追踪网络动态，使用EEG/MEG来描述那些依赖模块化的快速（亚毫秒级）波动。
本文假设在静息状态下，脑的动态模块化重组的特征是不同RSNs内部和之间进行分离和合并的连续过程。本文使用最近开发的框架来检测“模块化大脑状态”，该框架能精确量化随时间变化的“状态”波动。与其他聚类算法不同，该框架能检测大脑模块化结构的快速瞬时变化，并在时空精度方面优于其他现有聚类算法。在三个独立的静息态健康受试者的EEG/MEG数据集（N=568）上验证该假设，源重构（source-reconstructed）到整个皮层的68个区域。使用EEG/MEG源连接技术（利用功率耦合和相位耦合），结合滑窗方法和检测模块化状态的算法来重构动态脑网络（图1）。值得注意的是，本文结果表明了大多数受试者出现一致的网络模式，以及其中一些网络存在分裂现象，例如默认模式、视觉、颞区和背侧注意网络。本文推测，追踪正在进行的神经元活动的快速模块化架构，为人脑大规模电生理网络组织的动力学提供了新的见解。
图1 研究框架(A)分析三个数据集。(B)基于解剖图谱将模板MRI分割为感兴趣区域(ROI)。然后，对数据集1和2使用加权最小范数估计逆解（WMNE)和对数据集3使用波束形成来重建每个受试者的区域时间序列。(C)使用滑窗技术计算得到动态脑网络。(D)基于模块化算法的第一步是将每个时态网络划分为社区。然后，评估时间模块结构间的相似性。(E)相似性矩阵被分割成不同的社区，每个社区代表结合不同时间窗的特定空间拓扑的模块化状态。(F)提取每个受试者的所有模块（来自不同MSs）。只保留与RSNs呈强相关的（超过80%的节点重叠）。(G)计算与RSNs相关的主要模块的平均停留时间和分数占有率。
2. 结果
对三个独立的数据集进行分析：1)由健康大脑网络生物库提供的444名受试者的静息态脑电图数据；2)56名健康受试者静息态脑电数据；3)由人类连接组计划提供的61名受试者的静息状态MEG数据。对于EEG数据集（数据集1和2），使用加权最小范数估计(wMNE)，相位耦合。对于MEG数据集（数据集3），使用波束形成方法，包络耦合，校正空间泄漏以减少体积传导效应。然后，应用滑窗技术得到一系列连续的片段，表征每个人功能性脑网络的变化。随后应用基于模块化的算法，将动态网络的张量作为输入，并产生随时间波动的模块化状态(MSs)，每个MS代表一个特定的模块化拓扑。简言之，该算法通过量化所有时间模块间的相似性来确定具有相同拓扑的模块化结构。然后，确定与一个或多个RSNs强关联的模块。该研究的整体构架如图1所示。
2.1第一个数据集确定了16个状态
图2展示了从444名受试者的alpha波段中得到的16个模块，并报告了存在每个模块的受试者占比。得到三个与DMN相关的衍生模块：POST-DMN包括DMN的后部分成分（后扣带回、海马旁、楔前叶和顶下小叶区域），ANT-DMN包括DMN的前部分成分（前额叶区域）和DMN（将后部分和前部分整合到同一模块的大模块）。另外，得到三个颞区模块：L-TEMP和R-TEMP分别代表左右上颞区和下颞区；TEMP表示将左右颞区模块结合。还观察到与视觉网络相关的两个模块：ventral-VIS，包括视觉网络的腹侧区域；VIS整合腹侧和背侧视觉区域。总的来说，受试者百分比从高到低排列的模块是：DMN（存在于96%的受试者中）、POST-DMN（存在于88%的受试者中）、VIS-视觉网络（存在于86%的受试者中）、ANT-DMN（存在于86%的受试者）、SMN-躯体运动网络（存在于83%的受试者）、DAN-背侧注意力网络（存在于80%的受试者）、L-TEMP（存在于76%的受试者）、TEMP（存在于63%的受试者）、SAN-突显网络（存在于63%的受试者）、AUD+VIS-（结合听觉和视觉网络的一个模块，存在于61%的受试者中）、DMN+FPN-（结合默认模式和额顶网络的一个模块，存在于58%的受试者中）、DAN+VIS-（结合DAN和VIS的一个模块，存在于53%的受试者中）、R-TEMP（存在于45%的受试者）、FPN（存在于43%的受试者中）、AUD+VIS+DAN-（结合AUD、VIS和DAN网络的模块；存在于35%的受试者中）和腹侧-VIS（包括舌状和梭形视觉区域；存在于25%的受试者中）。结果在功能连接矩阵的几个阈值间和beta频带内是一致的。
图2提取数据集1中alpha波段的结果：与RSNs相关的衍生模块及其相应的受试者百分比。
2.2第二个数据集确定了12个状态
根据第二个数据集（图3），在alpha频段中，从57名受试者中得到了12个模块。这些模块是：POST-DMN（存在于98%的受试者中）、VIS（存在于94%的受试者中）、DAN（存在于91%的受试者中）、DMN（存在于84%的受试者中）、L-TEMP（存在于82%的受试者中）、ANT-DMN（存在于81%的受试者中）、SMN（存在于73%的受试者中）、AUD+VIS（存在于68%的受试者中）、DAN+VIS（存在于60%的受试者中）、TEMP（存在于45%的受试者中）、DMN+CCN（将DMN与认知控制相结合的模块，存在于32%的受试者中）和SAN（存在于22%的受试者中）。结果在功能连接矩阵的几个阈值间和在beta频带内也是一致的。
图3 提取数据集2中alpha波段的结果：与RSNs相关的衍生模块及其相应的受试者百分比。
2.3第三个数据集确定了10个状态
图4展示了数据集3中alpha波段的的结果，61名受试者确定了10个模块：DMN（存在于100%的受试者中）、POST-DMN（存在于95%的受试者中）、VIS（存在于88%的受试者中）、L-TEMP（存在于78%的受试者中）、SAN（存在于72%的受试者中）、DAN（存在于79%的受试者中）、SMN（存在于62%的受试者中）、AUD+VIS（存在于49%的受试者中）、DAN+VIS（存在于43%的受试者中）和CCN（存在于18%的受试者中）。结果在功能连接矩阵的几个阈值间和beta频带内是一致的。
图4 提取数据集3中alpha波段的结果：与RSNs相关的衍生模块及其相应的受试者百分比 总之，从三个数据集得到了与RSNs一致的波动模块。在所有数据集的受试者中，默认模式网络是最一致的网络（在受试者中表现出最高的百分比）。结果还表明，一些RSNs随时间推移而呈现出各种模块化拓扑结构，如DMN、颞区和视觉网络。此外，随着时间的推移而呈现出多个RSNs的模块组合，反映了跨网络交互。
图5显示了获得特定受试者结果的步骤，得到8个MSs。评估了从所有MSs和RSNs模板中提取的30个模块间的相似性（图5C）。在这30个模块中，有18个模块满足节点重叠大于80%的阈值条件。这18个模块与11个RSNs相关联：DMN、SAN、SMN、VIS、LTEMP、DAN+VIS、AUD+VIS、POST-DMN、DMN+FPN、FPN、ANT-DMN。通过观察图5C的相似性矩阵，显示DMN、SMN、VIS、ANT-DMN、AUD+VIS和LTEMP代表两个或多个模块，而SAN、DAN+VIS和POST-DMN与单个模块相关。图5D显示了所得到相关联模块的动态波动，其中每个模块根据其相应的RSN进行颜色编码。
图5特定受试者的结果。(A)提取特定受试者的8个模块状态的动态波动，(B)从每个MS中得到所有模块的空间表示。(C)所有MSs模块和RSNs模板间的相似性矩阵（仅提及与其中一个模块的重叠超过80%的RSNs），*标记重叠高于80%。(D)满足80%重叠阈值并与RSN模板相关的模块的动态波动。这些模块根据相应的RSNs进行颜色编码（如图5C右侧所示）。
2.4.停留时间和分数占有率
为了量化每个模块的时间特征，计算了两个指标：停留时间（DT），即在一个模块中存在连续窗的平均数量；和分数占有率（FO），即花费在每个模块中的时间比例。例如，图5中获得的模块按照FO排序是：DMN（FO=56%）、VIS（FO=43%）、DAN+VIS（FO=43%）、SMN（FO=31%）、SAN(FO=19%)、LTEMP(FO=14%)、AUD+VIS(FO=10%)、ANT-DMN(FO=4%)、POSTDMN(FO=3%)、DMN+FPN（FO=3%）、FPN（FO=3%）。根据DT，模块排序如下：DMN（DT=16%）、VIS（DT=10%）、DAN+VIS（DT=7%）、SMN（DT=7%）、AUD+VIS（DT=7%）、SAN（DT=6%）、POSTDMN（DT=6%）、LTEMP（DT=6%）、FPN（DT=5%）、DMN+FPN（DT=5%）。
图6展示了从每个数据集中获得的模块的FOs和平均DTs。在所有数据集中，DMN具有最高的FO和DT。从FO的角度出发，数据集2中的VIS网络较为重要。从DT的角度出发，数据集1中的SMN和SAN为非常稳定的模块。总之，从所有数据集获得的结果都证实了DMN的重要性（和稳定性）及在静息状态下作为功能性核心网络。
图6 小提琴图显示了分别从(A)Dataset1、(B)Dataset2和(C)Dataset3中得到的衍生模块的分数占有率（FO）和平均停留时间(DT)。每个图中出现的水平虚线表示平均值加两个标准差。∗标记重要模块（平均值>平均值+两个标准差）。
2.5.衍生模块与心理意象间的相关性
本文试图探究衍生模块与静息态获取期间所经历的心理意象（通过静息态问卷测量,rsQ）间是否存在相关性。只有数据集2有该数据。从rsQ中得到五个主要指标（即视觉心理意象、内心语言、体感意识、内心音乐体验和对数字的心理操纵）与每个模块的时间特征相关。任何rsQ分数和所提取模块的DT间不存在显著相关。
图7展示了视觉心理意象和从alpha波段得到的VIS、DAN和AUD+VIS的FO间呈正相关。在beta波段，结果显示视觉意象与AUD+VIS和DAN的FO之间呈正相关。
总之，这些结果表明，在采集过程中所经历的视觉意象的个体差异与特定模块的FO呈正相关，主要是VIS、DAN和VIS+AUD。
图7 视觉心理意象与在alpha波段中的VIS、AUD和DAN模块的FO显著相关（数据集2）
3. 讨论 本文展示了大规模电生理网络模块化结构的快速变化是如何塑造自发脑活动的。使用最近开发的算法来提取随时间变化的重复模块化脑状态。该框架应用于三个独立的EEG/MEG数据集，并证实了RSNs经历了连续模块化变化，这些变化存在于静息网络内部和之间的分离和合并的过程。 DMN动态切换其模块化拓扑，证实DMN实际上可分解为子成分，主要是前部和后部。研究还表明动态状态转换导致在一些大脑状态中包含FPN区域（数据集1的结果）。同样，颞区网络可分为左、右和完整模块。通过整合不同RSNs所构成的模块证实了静息脑的动态模块化行为。 本文研究结果表明，自发活动期间存在脑网络的动态变化。这些脑区是高度动态的，随时间的推移而改变其模块化排列。为了测试在大脑区域的一致性方面是否存在某种层次结构，本文报告了每个RSN所包含的区域及其在所有受试者中的占比（单位为%）。一旦模块与特定的RSN相关联，就会确定为重叠节点，以计算它们在所有受试者和数据集中的包含率。 DMN是在受试者中最一致的模块，因为它在受试者/数据集上的出现率最高，也具有高的FO值。说明DMN在整合自发脑活动中的信息方面起重要作用，因为存在大量与该网络相关的中枢。另外，DMN构型具有随着时间推移的最大稳定性（具有最长的平均停留时间）。停留时间短的模块是短暂的，在反映动态功能协调的时间内表现出大的可变性。本文研究中，衍生的短暂模块取决于每个受试者和数据集，大多属于高级认知网络、注意网络和感觉网络。在所有数据集中，这些停留时间较低的短暂模块是多个RSNs的整合（DAN+VIS、AUD+VIS…）。可能是因为在静息期间，个体动态地进行多种思维，而意象和走神仍是主要活动。因此，大脑是一个稳定活动与短暂功能变化相互交替的动态系统。 此外，本文研究强调了视觉网络的显著存在，显示出一段时间内的高占用率（数据集1和2的结果）。这种观察可能与大多数受试者在静息态采集期间表现出的视觉意象活动占主导地位有关。在采集过程中所经历的视觉意象的个体差异与VIS网络的FO呈正相关（图7）。心理意象与AUD、DAN和VIS之间的显著性关系可能说明这些网络间的交互作用，随时间的推移而形成一个大模块（在数据集1中观察到的结果，DAN+VIS+AUD的出现，图2）。 在三个数据集中，显示了9个共有的RSNs：DMN、POST-DMN、VIS、SMN、SAN、LTEMP、VIS+AUD、DAN+VIS。结果是基于记录信号的技术、预处理步骤、源重建、邻接矩阵阈值化、EEG/MEG频带（α波与β波）、图谱分割、功能连接测量以及校正或不校正零滞后相关性而获得。 其他模块分别来自每个数据集（图2、3、4）。每个衍生模块的受试者百分比也证实了受试者间的差异。在同一数据集中，这些个体差异与认知和行为功能的可变性有关。特定模块（主要是VIS、AUD和DAN）中的时间特征存在个体间差异与静息态问卷测量的心理视觉意象的自我报告评分有关。 虽然本文发现总体一致，但三个数据集的结果并不完全相同。如，DMN在数据集1、2和3上的存在率分别为96%、98%和100%，而VIS网络的存在率为86%、94%和88%。这些结果的差异可能与数据集本身存在的差异有关，如样本量、受试者年龄和实验条件（即闭眼/睁眼）。结果差异的另一个原因是使用不同的模态(MEG/EEG)。有研究证实，MEG/EEG的差异在研究瞬态静息态功能连接模式时较明显。 大量EEG/MEG研究证实了α至β频率范围（8-30Hz）在驱动自发大规模神经元相互作用中的相关性。本文研究得到以下几个结论：1）不同模块随时间波动，2）默认模式网络被认为是最一致、独立和稳定的模块，其模块化拓扑呈动态变化，3）在一段时间内观察到整合几个RSNs的模块反映了网络间的相互作用（DAN-VIS），4）一些特定模块的FO和心理意象呈正相关。 从方法论的角度出发，本文与先前处理相同数据集的研究使用相同方法（从数据处理到网络构建）。对EEG数据集使用wMNE/PLV组合来重建动态网络。对MEG数据集结合波束形成构造与带限功率包络之间的幅度相关（和正交化）。合适的窗宽是重建动态功能网络的关键问题。这取决于影响时间序列自由度的感兴趣频带和所使用的相关性度量。在EEG数据集中，在给定频段上选择提供6个“周期”的最小合适窗口长度。在MEG数据中，选择滑窗大小为0.5s的方法。为了提取快速瞬态模块，本文应用了基于模块化的算法来提取随时间推移的主要模块化脑状态。本研究使用的算法将大脑视为动态模块化网络来进行分割。4. 材料和方法4.1 EEG数据集4.1.1 数据集1（HBN）4.1.1.1 受试者 从健康大脑网络（HBN）数据库中，获取444名健康受试者的静息态脑电图数据。年龄在5-21岁之间。4.1.1.2 数据采集和预处理 在声音屏蔽室中使用128通道的EEG设备以500Hz的采样率采集高密度脑电数据，带通频率为0.1-100Hz。记录参考在Cz（头部顶点）。确保阻抗在40kOhm以下。采集5分钟的静息数据（闭眼状态）。使用Automagic Matlab工具箱对脑电数据进行预处理。多伪影抑制算法（MARA）用于检测和抑制伪影。选择每个受试者的4个40s无伪影时期。4.1.2 数据集24.1.2.1 受试者 共59名健康受试者（29名女性）。平均年龄为34.7岁（SD=9.1岁，范围18-55岁）。教育程度从10年学校教育到博士学位不等。没有服用任何药物，也没有神经或精神疾病。采集脑电数据之后，所有受试者都完成了静息状态问卷（ReSQ），由62个items组成，由五种主要类型的心理活动组成：视觉心理意象、内心语言、体感意识、内心音乐体验和数字的心里操纵。4.1.2.2 数据采集和预处理 每次脑电采集10分钟的静息数据（闭眼状态）。使用64通道的Biosemi ActiveTwo系统采集脑电数据采用标准的10-20导联系统。确保阻抗保持在20kOhm以下。预处理与之前研究一致。不良通道信号（完全平坦或被伪影污染）通过插值程序对半径5cm内的相邻电极进行恢复。保留电压波动在+80到-80的信号。选择每个受试者的4个40s无伪影时期。4.1.3 构建动态脑网络 对于两个EEG数据集，使用“EEG sourceconnectivity”方法结合滑窗方法来构建动态脑网络。“EEG sourceconnectivity”包括两个主要步骤：1）求解逆问题以估计皮层源并重构它们的时间动态，2）测量重构时序间的功能连接。 步骤如下： 1、通过使用Brainstorm确定解剖标志点来配准EEGs和MRI模板（ICBM152）。 2、使用OpenMEEG软件构建一个真实头模型。 3、基于Desikan-killiani图谱，将皮质表面分割成68个区域。 4、使用加权最小范数估计算法来估计得到区域的时序。 5、使用alpha波 8-13Hz和beta波 13-30Hz对重构区域的时序进行滤波。 6、使用相位锁定值（PLV）来度量重构区域时序间的功能连接。为了确保跨时间计算的所有动态网络密度相等，使用比例（基于密度）阈值方法来保留每个网络中前15%的连接值。4.2 MEG数据集（HCP）4.2.1 受试者 从HCP MEG2数据库中获得了61名健康受试者（38名女性）的静息态下的MEG数据。年龄从22-35岁。4.2.2 MEG采集和预处理 使用全脑Magens 3600扫描仪采集MEG数据。采集三个连续的session，每个session进行6分钟。要求受试者放松，睁着眼睛，并一直注视投射在黑暗背景上的十字准线。预处理包括去除伪影分割、确定错误记录的通道信号和回归掉在ICA分解中作为独立分量出现的伪影。4.2.3 构建动态脑网络 构建方法和以前处理同样数据集的方法一致。为解决逆问题，应用线性约束最小方差波束形成器。预先计算singel-shellsource models，分别计算1-30Hz和30-48Hz中的数据协方差。使用标准化导联场将数据波束形成到6mm网格上。然后，通过投影传感器噪声的功率对源估计进行标准化。源空间数据在alpha波1-13Hz和beta波13-30Hz中进行滤波。基于AAL模板得到区域时序后，使用对称正交化程序消除信号泄露。使用振幅包络相关测量（AEC）来估计区域时序间的功能连接。该方法简单概括为1）使用希尔伯特变化将功率包络计算为信号幅度，2）测量ROI功率包络的对数间的线性幅度相关性。最后，使用滑窗（length=6s，step=0.5s）构建动态连接矩阵。保留每个网络中15%的最强连接来进行阈值化。4.3 提取模块化脑状态 每个模块化状态反映了特定空间模块化组织。算法步骤如下： 1、运用Girvan-Newman和Louvain算法将每个时态网络分解成模块。每个模块算法重复200次迭代，对于每个网络产生400个模块化组织（200 runs2 algorithms）。运用Bassett提出的consensus algorithm来确定最终的模块化组织，该算法包括计算N的关联矩阵（N是节点数量），计算所得到的400个模块化组织中的两个节点被分配到同一模块的次数。然后，将关联矩阵与从原始分区的随机排列计算的零模型关联矩阵进行比较。仅保留关联矩阵的显著值。最后使用重复100次迭代的Louvain算法对阈值化关联矩阵进行聚类。 2、使用Rand系数的z-score（值介于0和1之间）来评估时间模块化结构间的相似性。这一步得到一个TT的矩相似阵，T是时间窗的数量。 3、使用consensus modularity方法对相似矩阵进行聚类，得到模块状态（MS）。4.4 提取与RNSs相关的模块 本文评估了每个模块（来自MSs）与RSNs间的相似程度。形成不同的masks或templates，每个mask与一个RSN、RSN的子部分或几个不同RSNs整合。然后，计算每个模块与每个模板的总体匹配。如果模块与模板间的重叠超过80%，就认为二者相关联。
4.5 量化 对与RSN模板相关的每个模块，计算了两个指标。 1、时间分数占比（FO），表示每个模块所出现的总时间，用百分比来衡量。因此，高FO值反映模块的高时间占比。 2、平均停留时间（DT）定义为在特定模块中所出现的连续窗口的平均数量。因此，具有高DT值被认为是稳定模块，低DT值被认为是瞬时模块。4.6 统计测试 通过评估静息状态问卷（rsQ）测量的认知表型间的统计关系，来研究所观察的大脑模块是否与静息态获得过程中的主观内在想法和感受有关。即计算rsQ中的五个主要指标（视觉心理意象、内心语言、体感知觉、内心音乐体验和数字的心里操纵）与衍生模块的FO之间的皮尔森相关。为了解决多重比较问题（5个心理活动与11个模块），使用Bonferroni校正p值，得到一个p<0.0009的调整后阈值。