导读

目的：人们对大规模脑网络功能连接动力学研究的兴趣不断增长，同时也在努力寻找电生理相关。在EEG数据分析中，应用于空域和谱域的常用约束可能会使部分神经活动无法识别。本研究提出了一种与BOLD功能网络连接动态相关的多模态EEG频谱模式的揭盲方法。

方法：利用平行因子分析方法对EEG频谱图进行盲源分离，现已被证明是揭示神经活动模式的一种有用技术。同时采集的BOLD fMRI数据采用独立成分分析法进行分解。利用滑动窗口相关法在成分的时间序列上计算动态功能连接，并根据相关系数值定义网络之间的连接状态。采用方差分析(ANOVA)检验评估网络间连接状态的动态与EEG频谱模式波动之间的关系。

结果：本研究发现了三种与网络间连接状态动态相关的模式。第一种模式在alpha、beta和gamma频段有主导峰，与听觉、感觉运动和注意网络之间的动态有关。第二种模式，在theta和低频alpha频段存在主导峰，与视觉和默认模式网络有关。第三种模式，同样在theta和低频alpha频段存在主导峰，与听觉和额叶网络有关。

意义：研究者在此前的研究中揭示了EEG频谱模式波动与大规模脑网络血流动力学之间的关系。在本研究中，研究者认为当EEG信号数据不受标准的空域和谱域约束时，这种关系在大规模脑网络之间的功能连接动态水平上也存在。

前言

人脑在功能上形成的大规模网络可以通过功能磁共振成像(fMRI)观察到。通过BOLD效应研究这些大规模脑网络(LSBNs)的最广泛的分析方法可能是通过独立成分分析(ICA)对数据进行盲源分离。功能连接(FC)是一种静态FC，并假设整个分析数据具有时间平稳性。在过去的十年里，人们对研究与神经元集群的同步和不同步相关的FC动力学越来越感兴趣。其假设FC只在短时间内稳定。计算动态功能连接(DFC)的一种常用方法是滑动窗口相关法，其关键约束是窗口大小。

人们一直在努力寻找静态LSBNs的电生理相关性；其他研究则探讨了动态FC与EEG功率波动之间的关系。处理这些同步数据集的电生理部分通常有许多限制。在空间信息方面，可以选择多电极平均信号或只选择电极子集信号。由于关注频谱特性，通常将EEG频谱划分为各种感兴趣的频段，而忽略了频率之间的关系。这种约束可能是LSBNs与EEG功率波动关系中存在被试间和被试内变异性的原因。

采用数据驱动技术对EEG数据进行多种盲源分离可以消除这一缺陷，从而发现EEG活动的多种模式。利用ICA对EEG频谱图进行两种分解，以空间谱域或时域模式表征LSBNs似乎是一种很有前途的方法。然而，它涉及到将一个三相数据(时×频×空)重塑为两种模式(如空间谱域×时域)，可以分解EEG频谱图中的任何三相依赖关系。

在对EEG频谱图进行多路分解的过程中，更进一步的方法是使用平行因子分析(PARAFAC)，它可以将一个三向EEG频谱图分解为一组三相(时×频×空)EEG活动模式。在之前的工作中，研究者证明了PARAFAC可以在EEG频谱图中发现和分离活动模式。有些模式具有神经生理学的起源，而有些则被报告为生理或技术伪影。在后续的工作中，研究者证明了PARAFAC分解的可靠性及其与LSBNs相关的时间稳定性。在本研究工作中，研究者证明了基于变分贝叶斯统计的PARAFAC也可以揭示与BOLD功能网络连接动态相关的时频空模式。揭示的几种模式在频谱特征的不同频段包含多个峰值。

方法

数据采集

本研究测量了50名健康被试(20名女性，30名男性，年龄19-40岁，平均年龄25岁)15分钟的静息态fMRI-EEG数据。要求被试在测量过程中闭上眼睛，躺着不动。正如被试在汇报过程中所报告的那样，没有人在扫描仪中睡着。该研究得到了当地伦理委员会的批准，所有被试在参与研究前都签署了知情同意书。

使用1.5T Siemens Symphony扫描仪采集MRI数据。获取t1加权MPRAGE结构像作为解剖参考(160个矢状面，矩阵大小为256×256，重采样至512×512，TR=1700ms，TE=3.96ms，FOV=246mm，FA=15°，层厚=1.17mm)。采用梯度回波测量静息态功能数据，回波平面成像300次(除小脑下部外全脑32个横断面，矩阵大小为64×64，TR=3s，TE=40ms，FOV=220mm，FA=90°，层厚=3.5mm)。使用32通道MR兼容的BrainProducts EEG系统在功能MR扫描过程中同时测量EEG数据。BrainAmp MR放大器与BrainCap MR(30个EEG电极，1个ECG电极和1个EOG电极)配以BrainAmp ExG MR放大器，可获得另外两个双极ECG信号和安装在被试胸部的加速度计信号。数据采样率为5kHz，并用BrainVision SyncBox同步MR采集和EEG采样。fMRI/EEG数据处理和分析管道如图1所示，并在后续段落中进行描述。

图1.fMRI/EEG数据处理和分析管道示意图。灰色/白色背景方块表示在被试/组水平上执行的步骤。利用独立成分时间序列上的滑动窗口相关性计算BOLD动态功能连接。将时变功能连接转化为网络间连接状态动态，并将其与3D(电极、时间、频率)EEG频谱图的PARAFAC分解所形成模式的时域特征进行比较。

EEG数据处理与分析

EEG数据的预处理和分析是在BrainVision Analyzer 2软件(BrainProducts，GmbH，Germany)中完成的，并在MATLAB 2014b(The MathWorks，Inc，Natick，USA)中运行。首先去除与MR相关的伪影。对于梯度伪影(GA)的抑制，本研究采用了一种基于模板减法的方法。模板是通过20 TR间隔的滑动平均值创建的。然后，在同步记录的GA过滤ECG数据中检测R波，并将其计时用于心冲击图(BCG)伪影的抑制。使用Butterworth 8阶零相位1-40Hz带通滤波器(48dB/八度)对数据进行滤波，降采样至250Hz，并重参考平均参考值。通过将每个电极的EEG数据分割为3s的周期(基于fMRI同步采集的TR)，利用快速傅立叶变换(FFT)将这些周期转换到谱域，并取傅立叶复系数平方的绝对值，创建PARAFAC模型的3D频谱图。因此，3D频谱图具有电极、时间和频率模式。

最后一个预处理步骤是每个电极和频率箱对应的所有时间序列的均值减去并缩放到单位标准差，以获得在PARAFAC模型估计中所有被试的这些模式的同等效果。将被试特定的3D频谱图沿时域模式连结，生成一组3D频谱图。使用PARAFAC模型对所有被试的3D频谱图进行分解。盲源分离基于变分贝叶斯统计，生成一组模式，每个模式具有空间、时间和频谱特征。根据研究者之前的工作和对这些数据进行PARAFAC模型估计的经验，研究者选择分解为10个模式。在10个PARAFAC模式(PRFC)中，丢弃了3个，因为它们揭示了同步fMRI-EEG记录中出现的各种技术伪影。在随后的分析中，使用了其余7个PRFC模式的时域特征来评估与BOLD功能网络连接动态的关系。

fMRI数据处理与分析

数据处理和分析在SPM 12(Wellcome Trust Centre for Neuroimaging，UCL，UK)、GIFT v4.0a(Mind Research Network，Albuquerque，USA)和MATLAB 2014b(The MathWorks，Natick，USA)环境中运行。预处理管道由一系列步骤组成，首先将所有功能扫描重新对准到第一次扫描，然后用RETROICOR滤除生理噪声，由此从同步记录的双极ECG和加速度计信号中提取所需的ECG R波和呼吸峰时间。由于TR=3s，还应用了时间层校正，空间归一化将数据转换为标准的立体定向MNI空间。

根据Barton等人(2015)和Weissenbacher等人(2009)的建议，回归出24个运动参数时间序列，5个从脑室提取的时间序列，4个从白质提取的时间序列，以抑制功能数据中残留的生理噪声。然后，使用Mask Explorer工具控制空间异常的数据。最后一个预处理步骤是在全宽半极大值(FWHM)中使用8mm高斯核卷积对功能扫描进行空间滤波。采用数据驱动组ICA来识别LSBNs。每个体素中的平均时间序列被缩放到100(GIFT中的强度归一化选项)。使用ICASSO程序进行10次迭代，以获得由InfoMax算法估计的最可靠的20个独立成分(IC)。然后使用GICA3技术对被试特定和网络特定的BOLD波动所对应的时间过程进行反向重建。

功能连接动态和网络间连接状态

在计算特定LSBNs之间的DFC之前，丢弃了5个明显表示伪影的ICs。然后在滑动矩形窗口中计算每个被试的15个IC时间进程的所有独特组合之间的Pearson相关系数，从而形成特定被试目的DFC时间序列。为了评估窗长对BOLD功能网络连接动态和EEG功率波动之间关系的影响，本研究计算了四种不同窗长(45、60、75和90s)的整体分析。

本研究定义了三种网络间连接状态(BNCS)：正向同步(PS)、反向同步(NS)和非同步(D)，它们对应于DFC时间序列中的高正、高负和接近零的相关系数。为了识别这些BNCS，对被试特定的DFC时间序列进行阈值化。根据所有相关系数(参见图4)分布中的局部极小值确定特定阈值，这些相关系数在所有被试的15个IC时间进程的所有组合之间形成了DFC时间序列。

网络间连接状态与EEG功率波动

为了获得BNCS时间序列和PRFC时域特征的相同时间，对PRFC时域特征执行了滑动窗口平滑。在相同长度的矩形窗口中计算均值，就像在DFC中计算相关系数一样。将PRFC平滑时域特征延迟(6s)到BNCS时间序列(通过丢弃BNCS时间序列的前两个值和PRFC平滑时域特征的最后两个值)，以考虑血流动力学响应延迟。

然后采用多因素方差分析来评估由PRFC平滑时域特征所代表的EEG功率波动与由BNCS时间序列所代表的网络间连接状态之间的关系。测试了跨被试连结的每个PRFC平滑时域特征均值的因素影响，其中第一个因素是跨被试连结的特定BNCS时间序列(固定效应)，第二个因素是被试者因素(随机效应)。BNCS时间序列和PRFC平滑时域特征的所有组合的多个多因素方差分析检验形成了一个105×7个p值的数组，通过错误发现率(FDR)在p＜0.05水平进行多重比较校正。对DFC计算的每个窗长分别执行FDR。

结果

EEG数据的PARAFAC分解

图2显示了PARAFAC分解得到的10种模式的空间和频谱特征。特定模式的3D频谱图数据可解释方差由同一图中的箱线图表示。所有10种模式可解释方差的中位数为5.2。

图2.基于50例健康被试EEG数据的3D频谱图PARAFAC分解结果。

分解的结果可以分为两组。第一组包括与MR环境中不可能完全抑制的伪影相关的模式。PRFC6的频谱特征包含多个峰值。PRFC3和PRFC8与MR梯度伪影的残留有关。这些PRFC的频谱特征在33Hz左右的伽马频段包含尖峰。对梯度伪影抑制前后的标准EEG功率谱进行视觉检查，发现残留的峰仍然存在，即32层TR 3s MR采集序列的高谐波为10.7Hz。第二组包含被认为起源于神经活动的模式。只有这些被用于后续分析(PRFC1，2，4，5，7，9和10)。最后一步对所有EEG数据进行预处理，PRFC模式的估计特征的正或负值可以解释为相对于平均功率基线的功率增加或减少。

fMRI数据的ICA分解

在fMRI数据分析部分，使用ICA将预处理的功能数据分解为20个可靠的独立成分。空间映射如图3所示。在目视检查了每个特定的成分(其空间映射和时间进程)之后，根据Griffanti等人(2017)研究中的建议发现了5个伪影成分。IC10的空间映射包含了一个典型的由运动引起的位于大脑边缘的环。IC1、2、5和12的空间映射位于脑室、血管和白质。在随后的分析中丢弃这5个成分。其余15个成分被认为是LSBNs的代表，根据Allen等人(2011)研究中的结果，通过视觉比较空间映射并对其进行功能特化。

图3.50例健康被试fMRI数据的组ICA分解结果。

网络间连接状态

基于从DFC计算中得到的所有相关系数的分布，定义了三种FC状态(NS、D和PS)(图4)。所有四个窗长的分布形状非常相似，其形状的局部极小值作为阈值(阈值1=-0.74，阈值2=0.80)。从这个意义上说，NS状态对应于由区间<-1.00，-0.74>的Pearson相关系数描述的两个LSBNs之间的功能连接；类似地，D状态的区间为(-0.74，0.80)，PS状态的区间为<0.80，1.00>。

图4.75s窗长的动态功能连接计算中所有相关系数的分布，并突出显示了网络间连接状态定义的阈值。阈值1=-0.74，阈值2=0.80。

LSBN连接状态及相关EEG功率波动

多元方差分析结果表明，三种PRFC模式的相对功率波动与一些LSBNs之间的BNCS动态有关。图5显示了经过FDR多重检验校正后具有统计学显著性的关系(p＜0.05)。由于关系的强度依赖于DFC计算的窗长，这里呈现的是在FDR之后至少有一个窗长显著的结果。

图5.PRFC模式的EEG功率波动与BOLD网络间连接状态之间具有统计学显著性的关系。

PRFC2的频谱特征在alpha和beta频段有优势功率增加，而在gamma频段优势功率减少。对于该模式的功率波动，结果发现了三种与BNCS动态的显著关系，其中包括听觉和注意网络(IC3，IC16)、感觉运动和注意网络(IC8，IC20)以及两个注意网络(IC15，IC20)。前两者在FDR后的75s和90s的DFC计算窗长变化显著；第三种关系仅在75s窗长中变化显著。事后方差分析检验(图6)显示，所有三种关系与PRFC2的操作方式相同。当alpha和beta频段的功率增加而gamma的功率减小时，上述LSBNs之间的功能连接从同步状态(无论是正向还是反向)变为去同步状态。

图6.方差分析事后检验显示了PRFC模式的相对功率波动均值，特别是BNCS，BNCS动态和EEG功率波动之间存在显著关系。

PRFC7频谱特征包含theta和低频alpha的主频率峰值。其功率波动与视觉和默认模式(IC6、IC14) BNCS动态密切相关。事后检验表明，这两个LSBNs的同步伴随着theta和低频alpha功率增加。经FDR校正发现，在窗长为75s的DFC计算中存在统计显著性。最后，FDR校正发现PRFC10频谱特征的功率波动与具有听觉和额叶功能特化的LSBNs的BNCS动态(IC3，IC18)之间存在显著关系。随着theta的增加和PRFC10频谱特征描述的低频alpha的减少，LSBNs的正向同步从去同步化变化为反向同步。

结论

人们对在电生理数据中寻找功能连接动态标志的兴趣日益浓厚，这表明EEG数据在处理过程中往往会受到空间和谱域的限制；这可能会削弱数据潜力。本研究提出了一种替代方案，即采用PARAFAC模型，一种不受这些约束或限制的方法。PARAFAC分解是一种有前途的工具，它能够揭示EEG频谱图中具有生理意义的活动模式，证明了分解在时间上和在被试间是稳定的，并揭示了与BOLD LSBNs相关的模式。总而言之，本研究引入了PARAFAC技术，该技术也能够发现反映BOLD功能网络连接动态的时频空模式。

原文：Spatial-temporal-spectral EEG patterns of BOLD functional network connectivity dynamics.

https://doi.org/10.1088/1741-2552/aab66b