EEG微状态预测并发fMRI动态功能连接状态

前言

静息态功能磁共振成像(rs-fMRI)测量的大脑功能连接在多个时间尺度上有所不同，并确定了循环的动态功能连接(dFC)状态。这些发现与不同的认知和病理状态有关，有可能作为疾病的生物标志物，但它们的神经基础仍然存在争议。在静息态EEG研究中也发现了一些重复出现的微状态，这些微状态被认为是反映时变大脑状态的大规模功能网络的准同步活动。在此，本研究假设fMRI衍生的dFC状态可能与这些EEG微状态有关。为了验证这一假设，本研究定量评估了EEG微状态预测从健康被试在休息时收集的同步EEG-fMRI数据中并发fMRI dFC状态的能力。通过训练随机森林分类器，本研究发现四种典型EEG微态预测fMRI dFC状态的准确率为90%，明显优于其他EEG特征，如频谱功率。本研究结果表明，EEG微状态分析产生了fMRI dFC状态的可靠特征，为dFC的电生理基础提供了证据，同时也进一步支持了EEG微态能够反映大规模脑网络的动态。

前言

静息态功能磁共振成像(rs-fMRI)测量的血氧水平依赖(BOLD)信号的自发波动已广泛用于研究大脑内在功能连接(FC)，揭示了许多所谓的静息态网络(RSNs)。有趣的是，最近发现FC在几秒到几分钟的短时间尺度上波动，通常被称为动态功能连接(dFC)，可能反映了大脑网络在响应内部和外部刺激时的不断重组。事实上，dFC已被证明与大脑状态(刺激/任务、闭眼、警觉、唤醒、睡眠、麻醉、药物操作)以及年龄、性别或疾病相关，并且它也被证明能够比静息态FC更好地预测大脑状态和疾病。一些研究进一步使用了模式识别技术，特别是聚类、主成分分析和字典学习，以识别反复出现的dFC模式，并可能反映出大脑潜在的状态。

尽管对dFC，尤其是dFC状态的研究产生了极大的兴趣，但它们的神经元基础仍然是一个悬而未决的问题。解决这个问题的一种可能方法是通过同时记录脑电图(EEG)和BOLD-fMRI来识别dFC的电生理相关性，从而利用EEG与神经元活动更直接的关系。有研究初步考察了dFC波动的EEG相关性，而非dFC状态：有研究首次聚焦于与不同唤醒状态相关的任务负默认模式网络(DMN)和任务正背侧注意网络(DAN)之间的连通性变化，并测量EEG功率在α和θ频段的波动。作者发现DMN-DAN连接降低反映了高唤醒期，这与之前的研究一致。Tagliazucchi等人(2012)也采用了类似的方法研究EEG功率与dFC波动之间的关系，发现清醒被试的特定脑区之间存在正相关和负相关，而在处于警觉状态之间的被试中，这些区域都有所变化。此外，dFC波动与根据图论计算的网络结构时间指数相关。这两项研究侧重于一些公认的EEG节律中的功率波动。此外，另外两项研究考察了超低频(＜0.1Hz)EEG波动及其与dFC的相关性，确定了休息期间功能连接的准周期时空模式。迄今为止，只有一项研究专门调查了dFC状态的电生理基础，该研究分析了睁眼和闭眼状态下的同步EEG-fMRI数据。在不同的dFC状态中发现了不同的EEG频谱特征，其中一种状态仅发生在闭眼状态，与EEGα功率的增加相关，进一步支持了警觉对dFC波动的影响。这些初步但有希望的结果促使我们探索更多的、信息量更大的EEG特征作为dFC状态的相关因素。

有趣的是，静息态EEG数据也显示出波动，这种波动可以用有限数量的头皮电位地形图来进行描述，这些电位在短时间内(~100ms)保持稳定：即所谓的微状态。在对健康被试的多项研究中，一致报告了四种典型的微状态，捕获了大量的EEG方差。与其动力学表现出的无标度特性相关联，发现EEG微状态与BOLD-fMRI RSNs相关。

在这里，本研究记录了健康被试在休息期间的同步EEG-fMRI数据，考察了EEG微状态预测并发fMRI检测到dFC状态的潜力，并将其与Allen等人(2018)报告的EEG频谱特征相比较。

方法

本研究由瑞士沃州人类研究伦理委员会(CERVD)批准，参与研究的9名健康被试(23±3岁，5男/4女)提供了书面知情同意书。该数据集是在瑞士洛桑EPFL生物医学成像中心获得的。

EEG-fMRI数据采集

使用8通道接收射频线圈，在主动屏蔽7T Magnetom扫描仪(Siemens，Erlangen)上进行成像。使用2D同步多层切片(SMS)梯度回波平面成像(EPI)序列(3×SMS和2×平面内GRAPPA加速)获得功能图像，TR/TE=1000/25ms，69个矢状切片和2.2mm各向同性体素大小，实现全脑覆盖。使用反向相位编码方向的短EPI采集(5 volumes)进行图像失真校正。采用T1加权3D梯度回波MP2RAGE序列获得全脑1mm各向同性结构图像。

使用两个与MR兼容的32通道BrainAmp MR plus放大器(Brain Products，Germany)记录EEG数据。采用BrainCap MR模型(EasyCap，Herrsching，Germany)，其中包含59个Ag/AgCl环形电极，按照10-20系统排列，一个专用电极用于参考，一个电极置于背部用于心电图(ECG)记录，四个改进的EEG电极用作运动传感器。采样率为5000Hz，与扫描仪的10MHz时钟同步。使用MRI系统生理监测装置的呼吸带，以50Hz的频率记录呼吸痕迹。

对于每个被试，在MR扫描仪内同时进行8分钟的EEG-fMRI扫描。被试在休息时睁开眼睛进行扫描，同时注视出现在MR兼容LCD屏幕(Cambridge Research Systems，Rochester，UK)中央的红色十字，通过放置在RF线圈上的镜子进行观察。这样，头部和眼球的运动被最小化。该装置利用了较高的BOLD灵敏度、更短的最佳TE和在7T时可用的高并行成像能力，产生了比现有的EEG-fMRI数据集(通常在3T时获得)更高时空分辨率的全脑fMRI数据。

MRI数据分析

预处理步骤

丢弃前10s的数据，使信号达到稳定状态，并使用FSL的工具BET去除非脑组织。随后，使用FSL工具MCFLIRT执行时间层和运动校正，使用反相编码采集后用FSL的TOPUP进行B0-unwarping步骤，以减少EPI失真。然后，采用截止周期为100s的高通时间滤波，并使用半值全宽(FWHM)为3mm的高斯核进行空间平滑。使用以下回归因子进行线性回归去除干扰波动(包括生理噪声)：(1)使用RETROICOR四阶傅里叶级数对与心脏和呼吸周期相关的准周期BOLD波动进行建模；(2)与心率以及呼吸深度和速率变化相关的非周期性BOLD波动通过与相应脉冲响应函数的卷积进行建模；(3)白质(WM)和脑脊液(CSF)掩膜中测得的平均BOLD波动；(4)MCFLIRT估计的六个运动参数；(5)扫描与头部大幅度运动期间获得的体积相关的零回归量；这些是使用FSL的程序FSL_motion_outliers确定的，因此首先计算Power等人(2012)研究中提出的DVARS指标，然后在第75个百分位数加上1.5倍四分位间距处进行阈值设置。

对于每名被试，使用ANTs工具Antropos将T1加权结构像分割为灰质、WM和CSF，从各自的T1加权图像中获得WM和CSF掩膜。使用FSL工具FLIRT将功能图像与各自的T1加权结构像进行共配准，随后使用FSL工具FNIRT将功能图像与蒙特利尔神经研究所(MNI)模板进行配准。将WM和CSF掩膜都转化为功能空间，然后使用3mm球形内核进行侵蚀，以最小化部分容积效应。此外，根据Chang和Glover(2009)中描述的基本原理，侵蚀的CSF掩膜与MNI空间的大脑室掩膜相交。

根据自动解剖标记(AAL)图谱，将每个被试的结构图像分割成R=90个大脑非重叠感兴趣区域(ROI)。这些ROI被共配准到被试者的功能空间，预处理后的BOLD数据在每个ROI内平均。所得信号以0.1Hz的截止频率进行低通滤波，因为神经元起源的同步BOLD波动主要发生在该频率范围内。

dFC状态的估计

根据Preti等人(2017)和Abreu等人(2019)的建议和结果，使用35.0s(35TRs)的窗长和5.0s(5TRs)的步长，通过滑动窗口相关法估dFC。计算每个滑动窗口的所有ROI平均BOLD信号对的成对Pearson相关系数。对于每个被试，s(s=1, ..., 9)，最终的dFC矩阵，Cs，通过提取每个相关矩阵的上三角部分，将其向量化并减去该被试的静息态FC，得到

，其中M=(R2-R)∕2和Ts=88表示窗口的数量。将所有被试的Cs矩阵连接到一组dFC矩阵

，共有T=9×Ts=792个窗口。

为了识别dFC状态，本研究采用了L1范数正则化字典学习(DL)方法，遵循Abreu等人(2019)提出的方法。简而言之，这可以表述为矩阵分解问题C=DA，其中

和

分别表示dFC状态和相关权重时程；k为dFC状态的个数。通过求解所给出的优化问题来估计：

使C的重构误差

最小化；‖·‖F表示矩阵的Frobenius范数。使用MATLAB®工具箱SPArse建模软件(SPAMS)中实现的算法执行D和A的估计。在L1范数正则化框架下，利用非负参数λ控制解的稀疏性。通过改变k(从5到10，单位步长)和λ(从1到0.1259的十个值，递减指数步长)和最小化贝叶斯信息准则(BIC)来确定最佳k和λ值，这会降低模型的复杂性，因此有利于对dFC状态进行更简洁的估计。BIC是通过5-5倍嵌套交叉验证程序计算的，以解决潜在的乐观、过拟合偏差估计。为了进一步验证所提出的DL方法，研究者将其与最常用的dFC状态估计方法(k-means聚类和PCA)的BIC值进行了比较。

EEG数据分析

预处理步骤

EEG数据预处理包括以下步骤：梯度伪影校正；插值不良通道；时域带通滤波(1-45Hz)；脉冲伪影校正；降采样至500Hz；运动伪影校正；以及基于ICA的去噪。

EEG特征提取

从预处理的EEG数据中提取频谱功率和微状态特征，用于随后的fMRI dFC状态分类。为此，首先使用与dFC估计相同的滑动窗口法分割每个被试的EEG数据：窗长为35.0s，步长为5.0s。这允许从EEG中获取动态(逐窗)概要特征，与fMRI dFC特征的相关时间尺度相匹配。然后从每个窗口提取以下EEG特征。

频谱功率。在FFT变换后，计算每个窗口的每个通道中特定频带的频谱功率，然后用总功率进行归一化，以说明潜在的被试间变异性。此外，与dFC分析一致的是，减去了静息态频谱功率(所有窗口的平均值)；通过这种方式，可以评估EEG频谱功率波动(而非频谱功率本身)与fMRI dFC状态之间的关系，这遵循了之前关于fMRI连通性的基于频谱功率的E/MEG相关研究的基本原理。为了与Allen等人(2018)的工作进行直接比较，本研究还在不减去静息态频谱功率的情况下测试了基于频谱功率的模型。并且考虑与已知EEG节律相关的频带：δ(1-4Hz)、θ(4-8Hz)、α(8-13Hz)、β(13-30Hz)和γ(＞30Hz)。然后获得每个EEG窗口归一化频谱功率的空间分布(地形图)。

微状态分析。采用两种不同的方法对每个时间窗的EEG数据进行微状态分析：凝聚层次聚类法(AAHC；EEGLAB插件可在https://www.thomaskoenig.ch/index.php/software/10-EEGLAB-plugin-manual上获取)和地形时频分解法(TTFD)。对于AAHC，首先计算EEG全局场功率(GFP)，并检测其局部极大值，因为它们代表最高电场强度的瞬间，因此也代表地形信噪比(SNR)。随机选择GFP局部极大值的子集，然后将相关的EEG拓扑图用作聚类输入。然后确定时域(微状态，MS)中最主要的拓扑结构。对于TTFD，首先使用Garbor小波对多通道EEG数据进行时频分解；频率范围为1-45Hz，步长为1Hz。对于所得小波系数拓扑的实部和虚部，获得了3D数据集(通道×频段×时间)。这些数据在时间上连接起来，然后与最后两个维度的关系进行矢量化，得到最终的2D数据集(通道×[频段+2×时间])，从中计算GFP并提取其局部极大值。随机选择其中的一个子集，并将相关EEG拓扑图用作聚类输入，得出该时频域(MSTF)中最主要的拓扑图。通过将每个小波系数的地形(即2D数据集的每一列)分配给最相似的MSTF，可以获得表征每个MSTF的时频谱图。

对于这两种微状态估计程序，根据先前的文献(定义了大量健康被试的四种标准微状态)，将聚类数量设置为4。相比之下，GFP局部极大值的数量在1000、2500、5000和10000之间变化，并计算由所得微状态解释的每个窗口的EEG方差。对于两种微状态估计程序(AAHC为73%，TTFD为77%)，使用2500个拓扑图进行聚类时，发现EEG窗口的中值方差最大。因此，最终的MS和MSTF集由4个聚类和2500个输入拓扑图获得，对于所有EEG窗口，忽略拓扑图的极性。

为了量化每个微状态i=1，…，4对相应EEG数据的描述程度，全局解释方差(GEV)计算如下：

其中，ti表示分配给微状态i和Ci,T（ti）的时间点，即微状态i与ti处EEG拓扑图T之间的空间相关系数。GEV值被归一化，即∑iGEV(i)=1。

dFC状态分类

基于EEG的模型

为了根据EEG预测dFC状态，生成了以下模型，包括EEG频谱功率和微状态特征的不同组合(如图1所示)。

图1.EEG特征提取。

①频谱功率：与五种EEG节律(SPδ、SPθ、SPα、SPβ和SPγ，用总功率归一化并减去各自的静息态频谱功率)相关的平均频谱功率矢量拓扑图；每个特征向量的维数等于EEG通道数(59)。还考虑了一个连接所有矢量化拓扑图的模型SPALL，在这种情况下，产生了维数为295的特征向量(5个EEG节律×59个EEG通道)。

②微状态：使用AAHC/TTFD获得的微状态的矢量化拓扑图，分别具有最高的(MS1/MSTF1)、两个最高的(MS2/MSTF2)、三个最高的(MS3/MSTF3)和四个最高的(MS4/MSTF4；所有微状态)GEV值。对于MSi和MSTFi(i＝1，…，4)，每个特征向量的维数等于EEG通道数的i倍(即59、118、177或236)。

③微状态和各自的频谱功率：通过TTFD获得的具有最高GEV值的微状态矢量化地形，与随时间变化的相关平均频谱功率(MSTF1+SP1)相连。在这种情况下，每个特征向量的维数等于通道数加上频率bins(即59+45=104)。

由于微状态的GEV值越高，微状态对给定dFC窗口的EEG数据越具有代表性，因此假设高GEV值的微状态包含更多与dFC状态分类相关的信息。

分类

根据对相应dFC矩阵贡献最大的dFC状态标记每个时间窗；定义dFC状态标签的结果集为响应变量Y，定义解释变量X为基于EEG的特征向量(图2和3)。分类程序的目标是确定满足Y=f(x)的函数f。这里，使用随机森林进行分类；这类分类器首先构建大量去相关树，使用bagging的概念随机选择训练样本和特征用于分类。然后，通过对来自所有树的结果进行多数表决来获得树集合的最终分类结果。在这项工作中，训练了50棵去相关树的随机森林。随着学习到的树的数量增加，通过视觉检查树集合的累积预测误差来经验地确定树的数量，在所有检查的情况下，都达到了约50棵树的平稳水平。

图2.用于预测dFC状态的模型特征向量。

图3.定义分类问题的响应变量(dFC状态标签)和解释变量(基于EEG的特征)。

使用基于leaveone-subject-out(LOSO)框架的交叉验证(CV)计算结果分类器的准确性(ACC)。因为本研究的数据集包括9名被试，所以进行了9折CV：对于第s折(即第s个被试)，属于所有其他被试的Y和X条目用于训练，而属于第s个被试者的条目用于测试，得出每折的ACC值。对所有被试重复此程序，通过对折叠处的ACC值取平均值来获得总体ACC。对于每个类(即dFC状态)，提取其他指标，即精确率、错误发现率、召回率和漏报率。考虑到使用所测试模型的过拟合风险，还使用了Hastie等人(2009a)计算的相对过拟合率(ROR)的估计值：

 其中

表示预测误差的估计值(1-ACC)，

 

 为训练误差，并且

 

为无信息率。计算每折的ROR，然后对各个折叠进行平均。由于在训练随机森林分类器时随机选择训练样本和特征，分类过程重复了100次；然后通过对重复次数的折叠平均ACC和ROR值进行平均，得到最终ACC和ROR。

控制测试

在构建基于MS和MSTF的模型时，根据每个窗口的GEV对微状态进行排序的方法与基于空间相似性在所有EEG窗口匹配微状态的更传统方法进行了比较。此过程在EEGLAB工具箱的微状态插件中实现，并在每个窗口中迭代排列微状态，以获得最大的共通性，直到达到收敛。然后使用这些基于MS和MSTF的新模型重复分类过程。使用AAHC和TTFD算法，从整个EEG数据中估计四种微状态，也可以考虑最传统的微状态分析。对于每个窗口，计算这些微状态的GEV并用于排序，这是区分每个窗口dFC状态的唯一因素。除了ROR测量之外，本研究通过添加高斯白噪声和操纵噪声模型的SNR来扰动产生最高ACC的模型；SNR在0.1，1和10之间变化，并对每个噪声级重复分类过程。

结果

dFC状态的估计

通过DL以及k-means和PCA确定的dFC状态的嵌套交叉验证获得BIC值。与Abreu等人(2019)的结果一致，基于DL的方法明显优于k-means和PCA。特别是，用k-means初始化的DL在折叠中平均表现出最低的BIC值。在这种情况下，dFC状态的最佳数目为k=10，最佳正则化参数为λ=0.7943。采用最优方法获得的10组dFC状态如图4所示。

图4.dFC状态的估计。

EEG特征对BOLD-fMRI dFC状态的预测

接下来，提取EEG频谱功率和微状态特征，以预测fMRI dFC状态。如图1所示，从被试的EEG中提取特征。出于显示目的，这些特征代表每个dFC状态(而不是每个窗口/被试)，并通过对具有相同dFC状态标签的所有窗口/被试的特征进行平均。关于微状态，仅显示GEV最高的状态(MS1和MSTF1)；此外，仅显示了产生最高ACC的频带的平均频谱功率地形。当使用AAHC和TTFD算法时，Michel和Koenig(2018)研究中报告的四种标准微状态都出现了。MSTF谱图的平均谱功率随时间的变化大致遵循指数衰减，与控制EEG谱的幂律一致。有趣的是，具有前后偶极子的微状态在10Hz时显示出清晰的峰值，通常与α活动相关，这与Koenig等人(2001)研究闭眼时观察到的觉醒-静息态EEG记录一致。

基于EEG的dFC状态分类

表1描述了测试模型的ACC和ROR分类结果。基于微状态的模型产生了最高的ACC值，其中MSTF4的性能最佳(ACC=89.6%，ROR=0.01)，其次是MS4(ACC=33.2%，ROR=0.07)。其余模型(SP和MSTF1+SP1)表现不佳，ACC值低于40%。同样，Allen等人(2018)提出的不减去静息态频谱功率的模型也产生了较低的ACC值，在10%到20%之间。

表1

最佳模型MSTF4以及相应的混淆矩阵如图5所示。在混淆矩阵中，对角线项对应于正确分类的观察值(即dFC出现)，而非对角线项对应于分类错误的dFC出现。除了dFC状态#9之外，不同类(即dFC状态)之间dFC出现次数的分布相当均衡。关于MSTF4模型，除#7和#8之外，所有dFC状态的精度和召回率始终高于80%(错误发现率和漏报率低于20%)，甚至有一半的状态高于90%。在状态#7和#8的情况下，标记为dFC状态#7的dFC事件始终被错误地分类为属于dFC状态#8。当观察每个dFC状态的代表性微状态时，按其GEV排序(从上到下)，dFC状态#7和#8分类不佳的原因显而易见，因为其代表性微态在地形和相关的GEV方面是相似的。尽管MS4的ACC(83.2%)比MSTF4(89.6%)小，但MS4在所有dFC状态下的准确率和召回率都达到60%以上，而MSTF4在对dFC状态#7进行分类时的值低至33.3%，除了微状态B(采用文献中公认的术语)始终表现出最低的GEV值外，各微状态的地形和GEV值在一般dFC状态中发生了显著变化。对于MSTF4模型，不同dFC状态的微状态可变性被量化为每个dFC窗口的每个微状态与其余窗口的各个微状态之间的空间相似性(SS)。在相同标签的dFC窗口中，不同微状态和窗口的平均SS为0.65，而不同标签的dFC窗口的平均SS为0.50；根据单样本t检验发现，这种差异具有统计学显著性(p＜0.01)。

图5.具有最高精度的EEG模型（MSTF4拓扑图）。

在构建基于MS-和MSTF-模型的每个窗口上，通过GEV对微状态进行排序的方法与更传统的通过空间相似性对其进行排序的方法进行比较时，性能明显变差，ACC/ROR的值为减少/增加。其中，最佳模型MSTF4的ACC值为71.2%，与按GEV排序的MSTF4模型相比，下降了18.4%。尽管如此，使用该替代分类程序获得的ACC明显高于机会水平。相比之下，当考虑按GEV排序的4个MS和MSTF微状态的固定集合时，分类精度急剧下降，分别为19.1%和23.6%。

通过添加高斯白噪声扰动的EEG模型，其ACC显著降低，特别是对于更具挑战性的场景(SNR=0.1)，在所有情况下，ACC值都低于35%。当SNR=10时，模型的性能没有受到实质性影响。

结论

本研究时第一个评估EEG微状态预测BOLD-fMRI中自发dFC状态能力的研究。具体而言，本研究定量评估了此类EEG特征在分类问题中预测同步fMRI dFC状态的能力。结果发现，基于EEG微状态的模型能够预测fMRI并发dFC状态，分类精度高。特别是，在时频域估计的准确性高达90%。此外，基于EEG微状态的模型明显优于基于EEG频谱特征的模型。总的来说，本研究结果进一步支持了dFC状态的电生理基础，并建立了EEG微状态与大规模脑网络动态密切相关的生物标志物。

原文：EEG Microstates Predict Concurrent fMRI Dynamic Functional Connectivity States.

https://doi.org/10.1007/s10548-020-00805-1

小伙伴们点个“在看”，加

(星标)关注茗创科技，将第一时间收到精彩内容推送哦～