导读
静息态fMRI数据的动态功能连接(dFC)分析通常通过计算滑动窗口相关(SWC)来执行，然后进行k均值聚类，以便将每个窗口分配到给定的状态。使用合成数据的研究表明，k均值性能高度依赖于滑动窗口参数和信噪比。此外，被试间的异质性来源可能会影响组级聚类的准确性，从而影响dFC状态时间属性的测量，如停留时间和分数占有率。这可能会导致关于组间差异的虚假结论(例如，将临床人群与健康对照组进行比较时)。因此，当应用于多被试队列并探索实现聚类性能最大化的方法时，量化k均值估计dFC状态时间属性的能力是否重要？在这里，本研究探索在聚类之前使用降维方法，以便将高维数据映射到低维空间，为后续的聚类步骤提供显著特征。在将k均值聚类应用于编码数据之前，研究者评估了使用深度自编码器进行降维的效果。将这种深度聚类方法与使用主成分分析(PCA)、均匀流形近似和投影(UMAP)以及使用L1或L2距离对原始特征空间应用k均值的降维方法进行了比较。研究者使用合成数据集(来自多个异质被试的数据)对聚类性能进行了定量评估。结果发现，深度聚类提供了最佳的性能，而其他方法通常不足以捕获dFC状态的时间属性。然后，本研究展示了每种方法在真实人类被试数据中的应用，并表明降维方法的选择对状态时间属性的组级测量有显著影响。
前言
功能磁共振成像(fMRI)数据的功能连接(FC)分析被用于表征和量化大脑活动的时空模式。在该分析中，将脑区或节点之间的血氧水平依赖(BOLD)信号的统计相似性(例如相关性)用作构建FC网络的边缘。这有助于使用图论来量化全脑动态。近年来，动态功能连接(dFC)被广泛应用于研究静息态下大脑活动的时变组织。一种常见的方法是计算滑动窗口相关(SWC)，从而产生一组FC矩阵，然后可以将这些FC矩阵聚类成一组重复出现的FC模式或状态，如图1所示。根据每个dFC状态的稳定性和可变性，使用诸如停留时间(状态改变前处于状态的平均持续时间)和分数占有率(在给定状态下所花时间的比例)等统计数据来量化大脑的时空动态。
图1.SWC框架概述，包括从BOLD fMRI数据构建FC矩阵，然后进行聚类以确定dFC状态。
因此，dFC分析的一个重要假设是能够将大脑动态分割成具有特定空间相关模式的状态，并且这些模式在比扫描长度更短的一段时间内保持静止(类似于EEG微状态)。这通常在实践中通过使用聚类方法(例如k均值)来实现，以识别哪些加窗FC矩阵属于有限数量状态中的哪一个，并根据相似性指标对矩阵进行分组。由于N个节点的FC矩阵维数为N(N−1)/2，并且分割可以有多达数百个节点，因此将FC矩阵聚类为dFC状态是一个高维(无监督)学习任务。然而，此步骤通常使用k均值聚类(尽管其他研究使用了谱聚类或分层方法)来执行。
一些研究评估了SWC检测大脑活动动态变化的能力，以及特定窗口参数对其有效性的影响。然而，很少有人对聚类步骤的性能进行定量评估。一些研究使用合成或替代数据验证了dFC聚类方法。这些研究一致表明，k均值性能高度依赖于滑动窗口形状和长度以及信噪比。Gonzalez-Castillo等人(2015)使用任务态fMRI来加强认知状态之间的切换，将每个任务block(休息、记忆、视频、数学)视为不同的dFC状态，因此提供了真实数据中的地面真值。该研究结果表明，当在被试水平聚类且分割中有大量区域(＞100)时，准确率很高，但对于较少的分割和更短的滑动窗长时，性能会下降。
K均值已被广泛用于队列研究中的dFC分析，以评估一系列神经系统疾病(包括帕金森病、精神分裂症、路易体痴呆和自闭症)以及健康认知和睡眠的分数占有率和停留时间。在比较不同组间大脑时空动态的应用中，被试间的异质性来源，如血流动力学响应函数(HRF)的形状和噪声水平，可能导致被试间的差异掩盖了潜在的dFC状态。很明显，聚类步骤中的误差将导致分数占有率和停留时间等属性的测量不准确，因此可能得出关于组间差异的虚假结论。
因此，当应用于多被试队列时，评估k均值准确量化时空dFC模式及其转换统计量的能力至关重要。在对高维数据进行聚类的其他应用中，通常在聚类之前应用降维方法，以便将数据映射到低维空间。这为后续的聚类步骤提供了显著特征，减少了不相关或噪声特征的影响。本研究评估了在聚类dFC状态之前使用降维方法的情况(图1)，并建议在对编码数据应用k均值聚类之前使用深度自编码器进行降维。本研究将这种深度聚类方法与使用主成分分析(PCA)、均匀流形近似和投影(UMAP)以及使用L1或L2距离对原始特征空间应用k均值的降维方法进行了比较。使用了合成数据集对聚类性能进行了广泛的定量评估，这些数据集来自多个异质性被试的数据。通过聚类精度、提取状态的FC矩阵与真实状态的FC矩阵之间的相似性、以及分数占有率和停留时间测量的误差来衡量性能。在合成数据中发现，深度聚类优于其他方法。本研究展示了每种方法在真实人类被试数据中的应用，并表明降维方法的选择对状态时间属性的组级测量有显著影响。
方法
数据
本研究使用了合成和真实的fMRI数据。一般来说，SWC分析所获得的数据通常包括每个被试数分钟的静息态fMRI数据。节点由预先确定的结构或功能分割方案定义，或通过在组水平上应用独立成分分析(ICA)生成特定研究地图。然后通过计算每个区域内的BOLD信号平均值，构造出N个节点的T个时间点的时间序列。
合成数据
使用SimTB(https://trendscenter.org/software/simtb/)合成数据，在时空分离模型下模拟一组N个节点的BOLD活动。在该模型中，时间序列数据是由神经事件序列与HRF的线性卷积构造的。生成了由底层dFC状态的时间过程控制的数据，其中在给定时间占用的状态决定了每个节点的活动对所有其他节点的影响。在每个时间点，特定于状态的神经事件以一定的概率发生(设置为默认值0.5)。当节点中发生特定于状态的事件时，这对功能上连接到该节点的所有节点中的事件振幅具有加法或减法效应，由该时间步长处的dFC状态定义。除了这些特定于状态的神经事件外，随机发生的独特事件被添加到每个节点的时间过程中，代表大脑活动中自发的特定于节点的波动。为了验证这些独特事件发生的振幅和概率是不同的，从而创建了具有不同噪声水平的合成数据。
本研究旨在评估dFC聚类在应用于包含多被试队列时的性能，因此研究者创建了具有模拟FC时间进程的被试数据。对于每个模拟数据，从隐马尔可夫模型(HMM)中对dFC状态的潜在时间过程进行采样。这确保了状态变化之间的可变间隔。不同被试之间的dFC差异可能是由噪声水平(即神经噪声和测量噪声)和HRF形状等属性的个体差异造成的。因此，独特神经事件的概率和振幅、BOLD信号中添加的高斯噪声的振幅以及HRF的参数在不同的模拟数据中是不同的。HRF参数是从SimTB定义的默认分布中采样。本研究还旨在确保聚类性能独立于状态FC矩阵和状态转移矩阵，因此研究者创建了多个数据集，每个数据集包含多个被试。状态FC矩阵和HMM的转换矩阵是随机生成的，以便对每个数据集来说都是唯一的，但在数据集中的被试之间共享。
每个数据集的重复时间(即抽样率)TR为2s，总持续时间为270 TR(9 min)。图2显示了随机生成的dFC状态集(a)，随机生成的转换矩阵(b)，以及相应状态的时间进程(c)的示例。注意，这里并没有改变SimTB模型的功能，只是将生成批量合成数据的过程自动化了。
图2.随机状态、转移矩阵和时间进程的示例。
该模型用于生成训练数据集，以优化降维方法的参数，以及用具有不同噪声水平、节点数、状态数、被试数量和HRF模型的验证数据集来评估一系列实验条件下的聚类性能。
人类数据
为了应用于真实世界的数据，研究者从人类连接组项目(https://www.humanconnectome.org)发布的HCP1200中获得静息态fMRI数据。已经应用于这些数据的处理步骤包括Glasser等人(2013)的空间预处理和Smith等人(2013)的时间预处理。空间预处理包括校正由梯度非线性引起的空间畸变，头动校正，B0畸变全局强度归一化的校正，以及2mm全宽半最大值(FWHM)平滑。时间预处理包括高通时间滤波(＞2000s FWHM)、伪影和运动相关时间进程的回归，然后是时间中心化和方差归一化。预处理后，使用FSL中的MELODIC进行空间ICA分析，从而获得组水平分割。然后将ICA空间地图集映射到每个被试的BOLD时间序列上，得到每个个体的节点时间序列。研究者使用了N=50的脑区分割，并选择了没有明显质量控制问题的数据。这些数据的TR=720ms，持续时间为1200 TR(14分24秒)。
滑动窗口相关
由T个时间点和N个脑区组成的多变量时间序列，通过SWC转换为一系列FC矩阵。窗口用于为所有节点选择时间序列的一小段。然后，窗口按给定的步长时间移动，以提取给定被试的整个时间序列中相同长度的重叠片段。在合成数据中，研究者检验了矩形和锥形(汉明和汉宁)窗形状，以及步长为1TR(2s)的30-60s范围内的窗长。通过从精度矩阵中估计协方差来测量每个窗口的FC，用L1范数进行正则化，其中正则化参数λL1使用交叉验证对每个被试进行估计。
聚类和降维
由于FC矩阵是对称的，因此提取了上三角区并进行了向量化。因此，在加窗和向量化之后，BOLD数据被转换为X×Y矩阵，其中Y表示被试数量乘以每个被试的窗口数量，X表示成对相关性(等于N(N-1)/2)。然后，将该数据的每列分配到一个团簇。在聚类之前引入一个降维步骤，将X×Y矩阵转换为d×Y矩阵，其中d是低维表征中的维数。研究者测试了原始k均值与降维步骤后应用k均值的性能(图1)。使用的降维步骤是PCA、UMAP和深度聚类(自动编码器和k均值)。所有方法均采用相同的k均值过程。通过给每个窗口随机分配状态标签，将结果与随机聚类进行了比较。
由于降维方法需要参数调优，这里生成了包含50名被试、5个状态、标准HRF和中等噪声的训练数据集，并使用长度为40s的矩形窗口进行处理。这些数据用于优化每种方法的参数，使用网格搜索参数值以最大化聚类精度。采用参数值的粗网格来防止过拟合。分别用15、25和50个节点生成单独的训练数据集，以针对这些情景的不同输入数据维数重新调整聚类方法。
K均值
采用dFC分析中常用的k均值聚类方法。选择局部方差最大值的示例FC窗口，并对这些窗口对应的FC矩阵应用128次k均值重复(最多1000次迭代)，每个都用k均值++算法初始化。给出每个数据点与其最近质心之间的最小平方和的质心集，然后用于初始化所有窗口的k均值聚类(最大10000次迭代)。除了默认的欧氏(L2)距离指标外，本研究还测试了曼哈顿(L1)距离指标，由于其高维数而常用于dFC分析。
主成分分析(PCA)
将PCA应用于所有FC矩阵，然后将第一个p主成分作为k均值聚类的特征，其中选择p以最大化使用合成训练数据的聚类精度。每个分割区域的结果p如表1所示。
表一
均匀流形近似和投影(UMAP)
UMAP是一种非线性降维技术，通过构建高维图表征，将数据投射到低维流形上，然后优化低维图，使其结构尽可能类似于高维图。高维图的结构是根据到最近邻m的距离确定的。更高的m会导致低维投影，更准确地捕捉数据的全局结构，而不是保留到邻居的局部距离。必须选择的其他参数是低维子空间中的维数u，以及低维表征中的点之间的最小距离v。使用UMAP将所有FC矩阵嵌入到一个低维子空间中，然后对嵌入数据进行k均值聚类。如前所述，选择u、v和m的值是为了使合成训练数据的聚类精度最大化(表1)。
深度聚类
深度学习为神经成像分析提供了强大的工具，包括结构MRI中解剖结构或病变的分割，基于任务态fMRI数据中的认知状态标注，或从功能连接网络中进行临床诊断。尽管这些有监督的深度学习应用需要大量的地面真值数据进行训练，但自动编码器可以用作无监督应用的降维方法。
自动编码器是一种人工神经网络，通过低维编码层将输入数据复制到输出。这个低维编码层形成的瓶颈效应迫使网络提取特征，从而通过解码层再现原始数据。在这种情况下，输入数据被用作训练目标，因此自动编码器可以用于无监督聚类应用中的降维。在这里，研究者使用自动编码器作为数据驱动的方法来确定组水平的特征空间，允许聚类应用于低维编码层提供的显著特征。这个框架称为深度聚类。
提出的深度聚类框架包括在所有FC窗口上训练一个自动编码器，然后对编码数据应用k均值聚类(图3)。本研究使用了一个具有三个编码层和一个对称解码器的全连接自动编码器。使用Adam优化器训练权重，以最小化输入和输出之间的均方误差(MSE)，训练了100个epoch，批尺寸为50。隐藏层采用整流线性单元(ReLU)激活函数，低维层和输出层采用线性激活函数。如前所述，选择每层中的单元数是为了使合成训练数据的聚类精度最大化(表1)。
图3.自编码器架构。
实验
合成数据聚类
研究者改变了模型的参数和预处理步骤，以评估以下每种情况下的聚类分析性能：
a)被试数量：10、50和100。
b)分割区域数量：15、25、50。
c)状态数：3、5、7。
d)血流动力学响应函数(HRF)：对SimTB中提供的两种HRF模型进行测试；标准HRF和Windkessel-Balloon模型。
e)噪声：低、中、高噪声数据集是通过改变基础神经时间进程中独特事件的概率和振幅而生成的，并将高斯噪声的振幅添加到BOLD信号中。
f)滑动窗口形状：矩形、汉明、汉宁。
g)滑动窗口长度：30s、40s、60s。
在每个数据集中，没有变化的参数设置为以下默认值：50个被试；25节点；标准HRF；高噪声；矩形滑动窗长为40s。为了验证聚类状态，通过配对对应FC矩阵之间的余弦相似性最大的状态，将聚类状态匹配为真状态。为了比较不同聚类方法之间的dFC状态质心，通过对属于该状态的所有FC窗口进行平均(而不是使用由k均值步骤导出的低维质心)来构造给定状态的代表性FC矩阵。聚类性能通过聚类精度(正确标记窗口的百分比)、聚类质心与真实状态之间的平均余弦相似性、停留时间(状态变化前占用的平均时间)和分数占有率(在给定状态下的总扫描持续时间的比例)的均方误差来衡量。为了计算停留时间MSE和分数占有率MSE，为每个被试的每个状态计算这些属性，然后计算这些测量值与真实值之间的平方误差，并在被试和状态中取平均值。对于每个参数集，研究者构建了5个数据集，并在每个数据集上执行每种方法的10次运行。然后，在给定的参数集上对每种方法的所有50次运行的性能指标进行平均。为了证明在每个数据集中跨被试的时间属性测量之间的差异，本研究从每种方法的一次运行中获得结果，并执行非配对双尾t检验，以检验与每个状态的分数占有率和停留时间测量的真实分布的显著差异。
聚类真实数据
从人类连接组项目数据集中选择了五组不重叠的100名被试。在每组中，应用矩形窗口长度为55 TR(39.6s)的SWC，步长为2 TR(1.44s)，每个被试有573个窗口。这是基于前人研究的基础上，结果表明在运动噪声不大的情况下，矩形窗口适合于检测dFC。在聚类示例FC窗口时，使用簇内到簇间距离的肘形判据为每个数据集选择簇数k。在所有100名被试的组中，确定了最佳聚类数为k=4。然后，使用表1所示的N=50个节点中相应的参数对每种方法进行聚类。为了进行后续的比较，通过对对应的FC矩阵之间的余弦相似性最大的状态进行匹配。
为了评估每种方法提供的状态时间属性测量值，研究者计算了每个被试中每个状态的占有率和停留时间。为了确定训练自编码器引入的随机性是否影响了重复运行深度聚类时状态时间属性的测量，对每个数据集执行了10次深度聚类分析，并使用单因素方差分析(ANOVA)来测试每个状态的分数占有率和停留时间运行之间的差异。为了确定降维方法的选择是否影响了状态时间属性的测量，对每种方法执行了一次运行，并使用单向方差分析来测试分数占有率和停留时间方法之间的差异。采用Benjamini-Hochberg方法进行错误发现率(FDR)校正。如果产生了显著的结果(FDR校正p＜0.05)，则使用未配对双尾t检验进行方法之间的事后两两比较。
结果
验证
为了检验每种降维方法的有效性，本研究使用合成数据量化了聚类性能。在检验的方法中，深度聚类在所有合成数据集中呈现出最高的精度(图4)。深度聚类方法在提取状态和真实状态之间具有最高的平均余弦相似性，最低的分数占有率MSE，以及在大多数数据集中的停留时间MSE最低(图4)。但在10个被试、15个节点和30s窗口的数据集上的结果除外，其中PCA和UMAP分别在平均余弦相似性和较低的停留时间均方误差上的结果略好。此外，在所有数据集中，深度聚类是在分数占有率和停留时间测量上优于概率的唯一方法。
图4.合成数据验证结果。
图5显示了应用于高噪声数据集的每种聚类算法结果的示例数据。可以看出，深度聚类的停留时间测量值的分布在统计上与地面真值没有区别。另一方面，对于使用L1距离的PCA或k均值的状态1、2和5，以及使用L2距离的k均值的状态1，这些结果与地面真值显著不同(图5c)。此外，使用L1或L2距离、PCA和UMAP从k均值测量的分数占有率与状态3的真值显著不同。状态3的高分数占有率测量结果表明，大量FC窗口被错误地分配给该聚类。这对状态3的状态FC矩阵有明显的影响(图5a)，其中大多数元素由于在大量FC窗口上平均而接近于零。相反，使用深度聚类得到的所有五种状态FC矩阵在视觉上具有与地面真值相似的结构。
图5.聚类分析结果来自应用于高噪声合成数据的每种降维方法的一次运行。
在真实数据中的应用
将每种方法应用于来自人类连接组项目中五组100名被试的真实数据。第一组的聚类结果如图6所示。除了UMAP之外，所有方法估计的状态FC矩阵在很大程度上都具有可比性，但与其他方法相比，深度聚类的状态4也存在差异(图6a)。状态1和状态3的分数占有率测量结果；状态2-4的停留时间测量结果存在显著差异(图6b；单因素方差分析，FDR校正p＜0.05)。图6c绘制了5个被试的状态时间进程，并展示了每种方法的差异。
图6.每种降维方法的聚类结果(来自人类连接组计划项目中的数据)。
结论
本研究证明深度聚类框架（包括用于在k均值聚类之前进行降维的自动编码器）可提高合成数据的dFC聚类性能。当应用于现实世界的数据时，与标准方法相比，这种性能的提高导致大脑状态的时间特征测量方面存在显著差异。这些差异定性地反映了在合成聚类结果中观察到的差异。
人类连接组项目数据：https://www.humanconnectome.org
SimTB模型：https://trendscenter.org/software/simtb/
本研究中使用的代码：https://github.com/apcspencer/ dFC_DimReduction
原文：Using deep clustering to improve fMRI dynamic functional connectivity analysis.
https://doi.org/10.1016/j.neuroimage.2022.119288
小伙伴们点个“在看”，加
(星标)关注茗创科技，将第一时间收到精彩内容推送哦～
静息态fMRI数据的动态功能连接(dFC)分析通常通过计算滑动窗口相关(SWC)来执行，然后进行k均值聚类，以便将每个窗口分配到给定的状态。使用合成数据的研究表明，k均值性能高度依赖于滑动窗口参数和信噪比。此外，被试间的异质性来源可能会影响组级聚类的准确性，从而影响dFC状态时间属性的测量，如停留时间和分数占有率。这可能会导致关于组间差异的虚假结论(例如，将临床人群与健康对照组进行比较时)。因此，当应用于多被试队列并探索实现聚类性能最大化的方法时，量化k均值估计dFC状态时间属性的能力是否重要？在这里，本研究探索在聚类之前使用降维方法，以便将高维数据映射到低维空间，为后续的聚类步骤提供显著特征。在将k均值聚类应用于编码数据之前，研究者评估了使用深度自编码器进行降维的效果。将这种深度聚类方法与使用主成分分析(PCA)、均匀流形近似和投影(UMAP)以及使用L1或L2距离对原始特征空间应用k均值的降维方法进行了比较。研究者使用合成数据集(来自多个异质被试的数据)对聚类性能进行了定量评估。结果发现，深度聚类提供了最佳的性能，而其他方法通常不足以捕获dFC状态的时间属性。然后，本研究展示了每种方法在真实人类被试数据中的应用，并表明降维方法的选择对状态时间属性的组级测量有显著影响。
前言
功能磁共振成像(fMRI)数据的功能连接(FC)分析被用于表征和量化大脑活动的时空模式。在该分析中，将脑区或节点之间的血氧水平依赖(BOLD)信号的统计相似性(例如相关性)用作构建FC网络的边缘。这有助于使用图论来量化全脑动态。近年来，动态功能连接(dFC)被广泛应用于研究静息态下大脑活动的时变组织。一种常见的方法是计算滑动窗口相关(SWC)，从而产生一组FC矩阵，然后可以将这些FC矩阵聚类成一组重复出现的FC模式或状态，如图1所示。根据每个dFC状态的稳定性和可变性，使用诸如停留时间(状态改变前处于状态的平均持续时间)和分数占有率(在给定状态下所花时间的比例)等统计数据来量化大脑的时空动态。
图1.SWC框架概述，包括从BOLD fMRI数据构建FC矩阵，然后进行聚类以确定dFC状态。
因此，dFC分析的一个重要假设是能够将大脑动态分割成具有特定空间相关模式的状态，并且这些模式在比扫描长度更短的一段时间内保持静止(类似于EEG微状态)。这通常在实践中通过使用聚类方法(例如k均值)来实现，以识别哪些加窗FC矩阵属于有限数量状态中的哪一个，并根据相似性指标对矩阵进行分组。由于N个节点的FC矩阵维数为N(N−1)/2，并且分割可以有多达数百个节点，因此将FC矩阵聚类为dFC状态是一个高维(无监督)学习任务。然而，此步骤通常使用k均值聚类(尽管其他研究使用了谱聚类或分层方法)来执行。
一些研究评估了SWC检测大脑活动动态变化的能力，以及特定窗口参数对其有效性的影响。然而，很少有人对聚类步骤的性能进行定量评估。一些研究使用合成或替代数据验证了dFC聚类方法。这些研究一致表明，k均值性能高度依赖于滑动窗口形状和长度以及信噪比。Gonzalez-Castillo等人(2015)使用任务态fMRI来加强认知状态之间的切换，将每个任务block(休息、记忆、视频、数学)视为不同的dFC状态，因此提供了真实数据中的地面真值。该研究结果表明，当在被试水平聚类且分割中有大量区域(＞100)时，准确率很高，但对于较少的分割和更短的滑动窗长时，性能会下降。
K均值已被广泛用于队列研究中的dFC分析，以评估一系列神经系统疾病(包括帕金森病、精神分裂症、路易体痴呆和自闭症)以及健康认知和睡眠的分数占有率和停留时间。在比较不同组间大脑时空动态的应用中，被试间的异质性来源，如血流动力学响应函数(HRF)的形状和噪声水平，可能导致被试间的差异掩盖了潜在的dFC状态。很明显，聚类步骤中的误差将导致分数占有率和停留时间等属性的测量不准确，因此可能得出关于组间差异的虚假结论。
因此，当应用于多被试队列时，评估k均值准确量化时空dFC模式及其转换统计量的能力至关重要。在对高维数据进行聚类的其他应用中，通常在聚类之前应用降维方法，以便将数据映射到低维空间。这为后续的聚类步骤提供了显著特征，减少了不相关或噪声特征的影响。本研究评估了在聚类dFC状态之前使用降维方法的情况(图1)，并建议在对编码数据应用k均值聚类之前使用深度自编码器进行降维。本研究将这种深度聚类方法与使用主成分分析(PCA)、均匀流形近似和投影(UMAP)以及使用L1或L2距离对原始特征空间应用k均值的降维方法进行了比较。使用了合成数据集对聚类性能进行了广泛的定量评估，这些数据集来自多个异质性被试的数据。通过聚类精度、提取状态的FC矩阵与真实状态的FC矩阵之间的相似性、以及分数占有率和停留时间测量的误差来衡量性能。在合成数据中发现，深度聚类优于其他方法。本研究展示了每种方法在真实人类被试数据中的应用，并表明降维方法的选择对状态时间属性的组级测量有显著影响。
方法
数据
本研究使用了合成和真实的fMRI数据。一般来说，SWC分析所获得的数据通常包括每个被试数分钟的静息态fMRI数据。节点由预先确定的结构或功能分割方案定义，或通过在组水平上应用独立成分分析(ICA)生成特定研究地图。然后通过计算每个区域内的BOLD信号平均值，构造出N个节点的T个时间点的时间序列。
合成数据
使用SimTB(https://trendscenter.org/software/simtb/)合成数据，在时空分离模型下模拟一组N个节点的BOLD活动。在该模型中，时间序列数据是由神经事件序列与HRF的线性卷积构造的。生成了由底层dFC状态的时间过程控制的数据，其中在给定时间占用的状态决定了每个节点的活动对所有其他节点的影响。在每个时间点，特定于状态的神经事件以一定的概率发生(设置为默认值0.5)。当节点中发生特定于状态的事件时，这对功能上连接到该节点的所有节点中的事件振幅具有加法或减法效应，由该时间步长处的dFC状态定义。除了这些特定于状态的神经事件外，随机发生的独特事件被添加到每个节点的时间过程中，代表大脑活动中自发的特定于节点的波动。为了验证这些独特事件发生的振幅和概率是不同的，从而创建了具有不同噪声水平的合成数据。
本研究旨在评估dFC聚类在应用于包含多被试队列时的性能，因此研究者创建了具有模拟FC时间进程的被试数据。对于每个模拟数据，从隐马尔可夫模型(HMM)中对dFC状态的潜在时间过程进行采样。这确保了状态变化之间的可变间隔。不同被试之间的dFC差异可能是由噪声水平(即神经噪声和测量噪声)和HRF形状等属性的个体差异造成的。因此，独特神经事件的概率和振幅、BOLD信号中添加的高斯噪声的振幅以及HRF的参数在不同的模拟数据中是不同的。HRF参数是从SimTB定义的默认分布中采样。本研究还旨在确保聚类性能独立于状态FC矩阵和状态转移矩阵，因此研究者创建了多个数据集，每个数据集包含多个被试。状态FC矩阵和HMM的转换矩阵是随机生成的，以便对每个数据集来说都是唯一的，但在数据集中的被试之间共享。
每个数据集的重复时间(即抽样率)TR为2s，总持续时间为270 TR(9 min)。图2显示了随机生成的dFC状态集(a)，随机生成的转换矩阵(b)，以及相应状态的时间进程(c)的示例。注意，这里并没有改变SimTB模型的功能，只是将生成批量合成数据的过程自动化了。
图2.随机状态、转移矩阵和时间进程的示例。
该模型用于生成训练数据集，以优化降维方法的参数，以及用具有不同噪声水平、节点数、状态数、被试数量和HRF模型的验证数据集来评估一系列实验条件下的聚类性能。
人类数据
为了应用于真实世界的数据，研究者从人类连接组项目(https://www.humanconnectome.org)发布的HCP1200中获得静息态fMRI数据。已经应用于这些数据的处理步骤包括Glasser等人(2013)的空间预处理和Smith等人(2013)的时间预处理。空间预处理包括校正由梯度非线性引起的空间畸变，头动校正，B0畸变全局强度归一化的校正，以及2mm全宽半最大值(FWHM)平滑。时间预处理包括高通时间滤波(＞2000s FWHM)、伪影和运动相关时间进程的回归，然后是时间中心化和方差归一化。预处理后，使用FSL中的MELODIC进行空间ICA分析，从而获得组水平分割。然后将ICA空间地图集映射到每个被试的BOLD时间序列上，得到每个个体的节点时间序列。研究者使用了N=50的脑区分割，并选择了没有明显质量控制问题的数据。这些数据的TR=720ms，持续时间为1200 TR(14分24秒)。
滑动窗口相关
由T个时间点和N个脑区组成的多变量时间序列，通过SWC转换为一系列FC矩阵。窗口用于为所有节点选择时间序列的一小段。然后，窗口按给定的步长时间移动，以提取给定被试的整个时间序列中相同长度的重叠片段。在合成数据中，研究者检验了矩形和锥形(汉明和汉宁)窗形状，以及步长为1TR(2s)的30-60s范围内的窗长。通过从精度矩阵中估计协方差来测量每个窗口的FC，用L1范数进行正则化，其中正则化参数λL1使用交叉验证对每个被试进行估计。
聚类和降维
由于FC矩阵是对称的，因此提取了上三角区并进行了向量化。因此，在加窗和向量化之后，BOLD数据被转换为X×Y矩阵，其中Y表示被试数量乘以每个被试的窗口数量，X表示成对相关性(等于N(N-1)/2)。然后，将该数据的每列分配到一个团簇。在聚类之前引入一个降维步骤，将X×Y矩阵转换为d×Y矩阵，其中d是低维表征中的维数。研究者测试了原始k均值与降维步骤后应用k均值的性能(图1)。使用的降维步骤是PCA、UMAP和深度聚类(自动编码器和k均值)。所有方法均采用相同的k均值过程。通过给每个窗口随机分配状态标签，将结果与随机聚类进行了比较。
由于降维方法需要参数调优，这里生成了包含50名被试、5个状态、标准HRF和中等噪声的训练数据集，并使用长度为40s的矩形窗口进行处理。这些数据用于优化每种方法的参数，使用网格搜索参数值以最大化聚类精度。采用参数值的粗网格来防止过拟合。分别用15、25和50个节点生成单独的训练数据集，以针对这些情景的不同输入数据维数重新调整聚类方法。
K均值
采用dFC分析中常用的k均值聚类方法。选择局部方差最大值的示例FC窗口，并对这些窗口对应的FC矩阵应用128次k均值重复(最多1000次迭代)，每个都用k均值++算法初始化。给出每个数据点与其最近质心之间的最小平方和的质心集，然后用于初始化所有窗口的k均值聚类(最大10000次迭代)。除了默认的欧氏(L2)距离指标外，本研究还测试了曼哈顿(L1)距离指标，由于其高维数而常用于dFC分析。
主成分分析(PCA)
将PCA应用于所有FC矩阵，然后将第一个p主成分作为k均值聚类的特征，其中选择p以最大化使用合成训练数据的聚类精度。每个分割区域的结果p如表1所示。
表一
均匀流形近似和投影(UMAP)
UMAP是一种非线性降维技术，通过构建高维图表征，将数据投射到低维流形上，然后优化低维图，使其结构尽可能类似于高维图。高维图的结构是根据到最近邻m的距离确定的。更高的m会导致低维投影，更准确地捕捉数据的全局结构，而不是保留到邻居的局部距离。必须选择的其他参数是低维子空间中的维数u，以及低维表征中的点之间的最小距离v。使用UMAP将所有FC矩阵嵌入到一个低维子空间中，然后对嵌入数据进行k均值聚类。如前所述，选择u、v和m的值是为了使合成训练数据的聚类精度最大化(表1)。
深度聚类
深度学习为神经成像分析提供了强大的工具，包括结构MRI中解剖结构或病变的分割，基于任务态fMRI数据中的认知状态标注，或从功能连接网络中进行临床诊断。尽管这些有监督的深度学习应用需要大量的地面真值数据进行训练，但自动编码器可以用作无监督应用的降维方法。
自动编码器是一种人工神经网络，通过低维编码层将输入数据复制到输出。这个低维编码层形成的瓶颈效应迫使网络提取特征，从而通过解码层再现原始数据。在这种情况下，输入数据被用作训练目标，因此自动编码器可以用于无监督聚类应用中的降维。在这里，研究者使用自动编码器作为数据驱动的方法来确定组水平的特征空间，允许聚类应用于低维编码层提供的显著特征。这个框架称为深度聚类。
提出的深度聚类框架包括在所有FC窗口上训练一个自动编码器，然后对编码数据应用k均值聚类(图3)。本研究使用了一个具有三个编码层和一个对称解码器的全连接自动编码器。使用Adam优化器训练权重，以最小化输入和输出之间的均方误差(MSE)，训练了100个epoch，批尺寸为50。隐藏层采用整流线性单元(ReLU)激活函数，低维层和输出层采用线性激活函数。如前所述，选择每层中的单元数是为了使合成训练数据的聚类精度最大化(表1)。
图3.自编码器架构。
实验
合成数据聚类
研究者改变了模型的参数和预处理步骤，以评估以下每种情况下的聚类分析性能：
a)被试数量：10、50和100。
b)分割区域数量：15、25、50。
c)状态数：3、5、7。
d)血流动力学响应函数(HRF)：对SimTB中提供的两种HRF模型进行测试；标准HRF和Windkessel-Balloon模型。
e)噪声：低、中、高噪声数据集是通过改变基础神经时间进程中独特事件的概率和振幅而生成的，并将高斯噪声的振幅添加到BOLD信号中。
f)滑动窗口形状：矩形、汉明、汉宁。
g)滑动窗口长度：30s、40s、60s。
在每个数据集中，没有变化的参数设置为以下默认值：50个被试；25节点；标准HRF；高噪声；矩形滑动窗长为40s。为了验证聚类状态，通过配对对应FC矩阵之间的余弦相似性最大的状态，将聚类状态匹配为真状态。为了比较不同聚类方法之间的dFC状态质心，通过对属于该状态的所有FC窗口进行平均(而不是使用由k均值步骤导出的低维质心)来构造给定状态的代表性FC矩阵。聚类性能通过聚类精度(正确标记窗口的百分比)、聚类质心与真实状态之间的平均余弦相似性、停留时间(状态变化前占用的平均时间)和分数占有率(在给定状态下的总扫描持续时间的比例)的均方误差来衡量。为了计算停留时间MSE和分数占有率MSE，为每个被试的每个状态计算这些属性，然后计算这些测量值与真实值之间的平方误差，并在被试和状态中取平均值。对于每个参数集，研究者构建了5个数据集，并在每个数据集上执行每种方法的10次运行。然后，在给定的参数集上对每种方法的所有50次运行的性能指标进行平均。为了证明在每个数据集中跨被试的时间属性测量之间的差异，本研究从每种方法的一次运行中获得结果，并执行非配对双尾t检验，以检验与每个状态的分数占有率和停留时间测量的真实分布的显著差异。
聚类真实数据
从人类连接组项目数据集中选择了五组不重叠的100名被试。在每组中，应用矩形窗口长度为55 TR(39.6s)的SWC，步长为2 TR(1.44s)，每个被试有573个窗口。这是基于前人研究的基础上，结果表明在运动噪声不大的情况下，矩形窗口适合于检测dFC。在聚类示例FC窗口时，使用簇内到簇间距离的肘形判据为每个数据集选择簇数k。在所有100名被试的组中，确定了最佳聚类数为k=4。然后，使用表1所示的N=50个节点中相应的参数对每种方法进行聚类。为了进行后续的比较，通过对对应的FC矩阵之间的余弦相似性最大的状态进行匹配。
为了评估每种方法提供的状态时间属性测量值，研究者计算了每个被试中每个状态的占有率和停留时间。为了确定训练自编码器引入的随机性是否影响了重复运行深度聚类时状态时间属性的测量，对每个数据集执行了10次深度聚类分析，并使用单因素方差分析(ANOVA)来测试每个状态的分数占有率和停留时间运行之间的差异。为了确定降维方法的选择是否影响了状态时间属性的测量，对每种方法执行了一次运行，并使用单向方差分析来测试分数占有率和停留时间方法之间的差异。采用Benjamini-Hochberg方法进行错误发现率(FDR)校正。如果产生了显著的结果(FDR校正p＜0.05)，则使用未配对双尾t检验进行方法之间的事后两两比较。
结果
验证
为了检验每种降维方法的有效性，本研究使用合成数据量化了聚类性能。在检验的方法中，深度聚类在所有合成数据集中呈现出最高的精度(图4)。深度聚类方法在提取状态和真实状态之间具有最高的平均余弦相似性，最低的分数占有率MSE，以及在大多数数据集中的停留时间MSE最低(图4)。但在10个被试、15个节点和30s窗口的数据集上的结果除外，其中PCA和UMAP分别在平均余弦相似性和较低的停留时间均方误差上的结果略好。此外，在所有数据集中，深度聚类是在分数占有率和停留时间测量上优于概率的唯一方法。
图4.合成数据验证结果。
图5显示了应用于高噪声数据集的每种聚类算法结果的示例数据。可以看出，深度聚类的停留时间测量值的分布在统计上与地面真值没有区别。另一方面，对于使用L1距离的PCA或k均值的状态1、2和5，以及使用L2距离的k均值的状态1，这些结果与地面真值显著不同(图5c)。此外，使用L1或L2距离、PCA和UMAP从k均值测量的分数占有率与状态3的真值显著不同。状态3的高分数占有率测量结果表明，大量FC窗口被错误地分配给该聚类。这对状态3的状态FC矩阵有明显的影响(图5a)，其中大多数元素由于在大量FC窗口上平均而接近于零。相反，使用深度聚类得到的所有五种状态FC矩阵在视觉上具有与地面真值相似的结构。
图5.聚类分析结果来自应用于高噪声合成数据的每种降维方法的一次运行。
在真实数据中的应用
将每种方法应用于来自人类连接组项目中五组100名被试的真实数据。第一组的聚类结果如图6所示。除了UMAP之外，所有方法估计的状态FC矩阵在很大程度上都具有可比性，但与其他方法相比，深度聚类的状态4也存在差异(图6a)。状态1和状态3的分数占有率测量结果；状态2-4的停留时间测量结果存在显著差异(图6b；单因素方差分析，FDR校正p＜0.05)。图6c绘制了5个被试的状态时间进程，并展示了每种方法的差异。
图6.每种降维方法的聚类结果(来自人类连接组计划项目中的数据)。
结论
本研究证明深度聚类框架（包括用于在k均值聚类之前进行降维的自动编码器）可提高合成数据的dFC聚类性能。当应用于现实世界的数据时，与标准方法相比，这种性能的提高导致大脑状态的时间特征测量方面存在显著差异。这些差异定性地反映了在合成聚类结果中观察到的差异。
人类连接组项目数据：https://www.humanconnectome.org
SimTB模型：https://trendscenter.org/software/simtb/
本研究中使用的代码：https://github.com/apcspencer/ dFC_DimReduction
原文：Using deep clustering to improve fMRI dynamic functional connectivity analysis.
https://doi.org/10.1016/j.neuroimage.2022.119288
小伙伴们点个“在看”，加
(星标)关注茗创科技，将第一时间收到精彩内容推送哦～
