《Frontiers in Computational Neuroscience》：Dynamic mode decomposition of resting-state fMRI revealing abnormal brain region features in schizophrenia
1 引言
精神分裂症是一种与大脑结构和功能异常相关的严重精神障碍。功能磁共振成像（fMRI）凭借其高空间分辨率，已成为研究这些神经改变的关键技术，因其提供的血氧水平依赖（BOLD）信号反映了局部脑能量消耗。目前，精神分裂症的诊断主要依赖临床访谈和评定量表，该方法因症状与其他疾病重叠、诊断标准模糊及解释主观性而面临挑战。因此，迫切需要客观、标准化的辅助诊断工具，神经影像学在此方面潜力巨大。
近年来，数据驱动算法（如多种机器学习模型）为神经科学提供了强大潜力，能够分析复杂神经系统并揭示超越传统方法的见解。其中，动态模式分解（DMD）是一种数据驱动方法，可从高维时间序列数据中提取具有物理意义的低维时空模式。其适用于非线性、高维、动态和时空多尺度及稀疏系统（这些特性内在于神经数据）的特性，使其非常适用于神经科学。除了理论上的适用性，DMD已被证明能有效从多种数据集中提取相干模式，包括神经生理记录。当应用于神经影像数据时，DMD提取的时空相干模式不仅补充了静态分析的发现，提高了神经解码的准确性，而且比时空独立成分分析（ICA）更能有效预测个体行为差异。这些优势确立了DMD作为表征脑活动时空组织的有前途的工具。
本研究将DMD应用于静息态fMRI数据，以提取全频段和分频段动态模式的平均振幅作为脑区特征。随后应用最小绝对收缩和选择算子（LASSO）回归选择关键特征，并利用支持向量机（SVM）验证其在精神分裂症患者与健康对照之间的分类效能。结果表明，通过多次LASSO选择获得的DMD特征能有效区分两组，并准确定位相应脑区，从而识别患者的异常脑区及振幅偏差。DMD不仅揭示了脑信号的动态演化模式，支持信号重建和预测，还证明了提取的特征与区域脑激活水平的相关性。这些发现为探索精神分裂症异常脑功能机制及识别精确诊断和治疗靶点提供了有效工具。
2 材料与方法
2.1 数据集与预处理
为研究DMD在区分精神分裂症患者与健康受试者的异常脑区特征提取和分类中的效果，本研究使用了生物医学研究卓越中心（COBRE）数据集。该数据集提供了72名精神分裂症患者和75名健康受试者（每组年龄范围18-65岁）的原始解剖和功能MRI数据。在去除两名脱落者和一名时间点缺失的受试者后，从剩余的145名受试者中选取一个子集，确保健康组与患者组人数相等，且在年龄、性别或利手方面无统计学显著差异（p > 0.05）。最终，分析纳入了N = 136名受试者的数据，包括68名精神分裂症患者和68名健康对照。
静息态fMRI数据使用DPABI-DPARSFA（V5.4_230110）进行预处理。预处理采用标准流程，包括去除前10个时间点、时间层校正、头动校正、去除头皮结构、将结构像对齐到功能像、分割、将图像匹配到欧洲模板、去噪、去除线性漂移、调整头动参数、滤波、归一化和图像平滑，未进行全局信号回归。使用‘Schaefer2018_300Parcels_7Networks_order_FSLMNI152_2mm’模板定义300个感兴趣区域，并提取每个受试者300个感兴趣区域的BOLD信号时间序列，计算脑区间的功能连接。
2.2 使用DMD进行特征提取
2.2.1 DMD算法概述
DMD算法的核心过程旨在从高维时间序列数据中提取关键动态模式。首先，对数据矩阵X进行奇异值分解（SVD）以实现降维和去噪，从而识别脑活动中最突出的空间模式。这为复杂数据建立了一个高效的“主坐标系”。随后，在这个简化空间中定义一个动态算子?，以捕捉系统状态随时间演化的核心规律。然后对?进行特征分解，解耦能够独立、稳定振荡的基本动态成分，每个成分由其增长/衰减特性和振荡频率表征。最后，将简化空间中识别出的这些动态成分重新映射回原始脑区空间，产生具有明确时空结构的DMD模式。每个模式对应一个完整的动态单元，在整个大脑中跨特定脑区表现出协调活动，并以独特的节律振荡。这个系统化的工作流程实现了从数据简化、动态原理提取到生成可解释神经特征的整个过程。
算法步骤包括构建数据矩阵X和X'（X在时间上向前移动一步），并假设系统状态的变换由矩阵A决定，即X' = AX。后续步骤包括对X进行SVD，计算简化空间中的动态算子?，对?进行特征分解以获取特征值和特征向量，最后将特征向量投影回原始空间得到DMD模式。
关键公式包括定义每个模式φi的模式振幅Pi为其向量模的平方（公式1），模式φi的时间振荡频率fi（公式2），以及使用初始状态x(0)、分解得到的模式φ（特征向量）和特征值Λ来重建时间t的BOLD信号（公式3）。公式3给出了静息态BOLD信号随时间变化的显式表达式。
2.2.2 数据矩阵构建与频率滤波
每个受试者预处理后得到的感兴趣区域时间序列矩阵的维度为n × m（n = 140个时间点，m = 300个脑区）。由于DMD算法的输入矩阵X应为特征×时间的形式，因此对每个受试者的时间序列矩阵进行转置后应用算法。每个受试者被分解得到139个模式，每个模式是一个长度为m（脑区数量）的复向量。
观察到数据中对应于0.01–0.1 Hz频带外频率的DMD模式振幅显著较小。这一经验观察与既定的神经影像共识一致，即BOLD信号在0.01–0.1 Hz范围内的波动（称为超慢活动）对静息态fMRI最具生理相关性，因为它们与自发神经活动和功能网络协调稳健相关。因此，该典型频带之外的模式被视为噪声并被消除。具体而言，高于0.1 Hz的频率主要受心脏（~1 Hz）和呼吸（~0.2–0.3 Hz）周期产生的非神经生理噪声污染，而低于0.01 Hz的超低频波动可能受扫描仪漂移或超慢生理过程干扰。应用此标准频率滤波器可增强解释内在脑动态的信噪比。
在0.01–0.1 Hz频带内，每个受试者有M = 62或64个模式，复模式以共轭对形式出现，每个模式有相应的频率。
2.2.3 特征提取与计算
Pi表征了该模式φi对应活动空间模式在整体信号中的强度，然后Pij（公式4）反映了相应模式φi在编号为j的各自脑区中的信号强度。由于模式是复向量，使用模式振幅是避免处理复数实部和虚部的好方法。计算每个模式i和每个脑区j的Pij，然后将每个脑区所有模式的Pij按模式数M平均，得到编号为j的每个脑区对应的平均信号强度featurej（公式5），将其作为每个受试者由DMD提取的全频段脑区平均振幅特征。
此外，为探索并利用DMD算法可同时提取时空信息的特性，还提取了每个受试子的子频段特征。即将0.01–0.1 Hz划分为若干频带，然后对每个频带内所有模式分量的Pij按该频带模式数平均，并将所有频带的平均值拼接在一起，得到每个脑区在不同频带的平均信号强度，作为每个受试者中DMD提取的子频段特征，用于探索精神分裂症患者脑区异常信号与频率的关系。
2.3 基于LASSO的特征选择结合SVM分类
为评估提取特征对受试者（患者/健康）的分类能力，采用SVM进行分类建模。通过比较不同特征集的分类性能，验证特征的有效性。在将特征输入SVM之前，对从每个受试者获得的脑区平均振幅特征（全频段和子频段）进行LASSO回归，以排除不重要的特征变量。LASSO回归的定义如公式6所示，其中y是目标变量，一个N维向量；X是大小为N × p的特征矩阵，其中N是样本数，p是特征数；β是大小为p维的回归系数向量；γ是正则化参数，控制L1正则化的强度。LASSO回归使用10折交叉验证进行拟合，初始化参数γ = 0.2。对于每一折，仅对训练集进行通过LASSO的特征选择。所选特征随后在相应的测试集上进行测试。通过计算准确率、精确率、召回率、F1分数和特异性这五个常用指标来估计分类器的有效性。在分类过程中，对数据集进行三轮十折交叉验证，产生30组性能指标。然后对不同特征的分类指标进行1000次置换检验，以判断它们是否具有统计学显著差异。
为验证通过LASSO选择的特征在区分患者组和健康组方面的功效，比较了三组脑区平均振幅特征的SVM分类性能：使用随机扰动训练标签获得的特征、未经LASSO选择的特征以及经过单轮LASSO选择后的特征。此外，由于单次LASSO回归具有固有的随机性，为实现稳定的特征选择结果，还比较了在三次独立的十折交叉验证过程中（即30次LASSO回归）通过1000次LASSO回归选择的交集脑区特征的SVM分类性能。这四组特征在相同的训练/测试集划分下进行三折交叉验证。通过比较它们的分类性能，全面评估了DMD提取的平均振幅特征和特征选择方法的鲁棒性和判别力。
考虑到单次LASSO回归可能因算法的固有随机性而产生不稳定结果（理论研究指出存在过早选择伪变量的风险以及倾向于过度选择特征的倾向），为获得绝对稳健的特征子集，进行了1000次独立的LASSO回归。最终特征集严格限制为仅包含在全部1000次运行中一致被选择的特征。此筛选标准克服了高维数据中变量选择的不稳定性，从而确保用于分类的最终特征集具有高重现性和统计可靠性。
2.4 精神分裂症患者异常脑区分析
将全频段和子频段脑区特征的平均振幅映射到相应的脑区，因为每个特征构成一个300维向量，对应300个不同的脑区。为充分利用现有数据探索精神分裂症患者的异常脑区，对所有受试者脑区平均振幅特征应用LASSO回归，识别对区分患者与健康对照至关重要的关键特征。然而，考虑到每次LASSO回归中变量保留的随机性，对全频段和子频段平均振幅特征进行了1000次独立的LASSO回归分析，以识别一致被选择的脑区特征（即无论进行多少次回归，在分类中均具有稳健主导作用的特征）。然后通过评估这些分析中的交集来确定一致保留的脑区。随后，对健康组和患者组之间的交集脑区特征进行了1000次置换检验。在显著性水平α = 0.05下显示出显著差异的脑区被定义为精神分裂症患者的异常脑区特征。
还计算了这1000次LASSO回归中这些异常脑区特征的平均回归系数β?。LASSO回归系数可以反映每个特征对预测结果的贡献和影响。具体而言，正值（β?i > 0）表示该特征与目标变量（患者标签为1，健康为0）呈正相关，负值（β?i < 0）表示负相关，绝对值β?i的大小表示影响程度。通过分析每个异常脑区的平均回归系数，可以了解精神分裂症患者异常脑区偏差的方向及其严重程度。
3 实验结果
3.1 使用全频段脑区平均模式振幅作为特征并执行LASSO回归的有效性
比较了四组不同的脑区平均振幅特征的SVM分类性能：使用随机扰动训练标签得到的特征（DMDp-chance）、未经LASSO选择的特征（DMDp）、通过单轮LASSO回归选择的特征（DMDp-LASSO）以及在三次独立十折交叉验证过程中通过LASSO回归识别的交集脑区特征（DMDp-LASSO-intersect）。
结果表明，在显著性水平α = 0.05下，DMDp、DMDp-LASSO和DMDp-LASSO-intersect的分类效果在所有指标上均显著优于DMDp-chance，且DMDp-LASSO-intersect的分类效果在准确率、召回率、F1分数和特异性上显著优于DMDp。DMDp-LASSO-intersect的分类效果在准确率和召回率上显著优于DMDp-LASSO。因此，使用脑区平均振幅作为分类特征，LASSO进行特征选择，以及将多次LASSO回归的交集作为关键判别因子，被证明是区分两组的有效策略。
3.2 精神分裂症患者的异常脑区活动
在全频段脑区平均振幅特征的1000次LASSO回归分析中，每次均保留了28个脑区。在显著性水平α = 0.05下，这28个脑区特征中有22个在健康组和患者组之间存在显著差异。将这22个具有显著差异的脑区特征视为在分类中起重要作用并能将精神分裂症患者与健康人群区分开的异常脑区特征。
绘制了这22个脑区（涉及脑区及编号如下：98-'顶叶’，108-'楔前叶’，249-'顶叶’，250-'顶叶’，262-'外侧前额叶皮层’，77-'额叶岛盖岛叶’，79-'额叶岛盖岛叶’，69-'中央前回腹侧’，210-'后部’，4-'视觉’，12-'视觉’，22-'视觉’，188-'躯体运动’，89-'颞极’，90-'颞极’，239-'眶额皮层’，244-'颞极’，246-'颞极’，123-'顶叶’，271-'顶叶’，275-'顶叶’，293-'背内侧前额叶皮层’）的平均回归系数β?的条形图以及以功能脑网络为分区的带误差棒的折线图。在SVM分类中，将患者标记为1，健康受试者标记为0。因此，LASSO回归的正系数表明患者组在该脑区的平均振幅高于健康组，负系数表明患者组在该脑区的平均振幅低于健康组，绝对值大小表示该脑区相对于健康组的异常程度。可视化展示了健康组和患者组在这22个异常脑区的平均振幅大小。
与顶叶、外侧前额叶皮层、额叶岛盖岛叶、中央前回腹侧、颞极、眶额皮层和背内侧前额叶皮层相关的脑区平均振幅在患者组中高于健康组；而与楔前叶、后部、视觉、躯体运动和顶叶相关的脑区平均振幅在患者组中低于健康组。属于显著性/腹侧注意网络和边缘网络的异常脑区在患者组中的平均振幅高于健康组；属于视觉网络和躯体运动网络的异常脑区在患者组中的平均振幅低于健康组。根据异常程度排名，前十个异常脑区分布在显著性/腹侧注意网络、视觉网络、躯体运动网络、边缘网络和默认模式网络中。
3.3 使用子频段平均振幅特征的分类效果及异常脑区
将0.01–0.1 Hz频带平均分为k个频带（k∈[2,8]），并将不同频带的平均振幅特征水平拼接，使得每个受试者有m?k个特征可用。在显著性水平α = 0.05下，将特征分为两个或三个频带所获得的分类结果显著优于未分割的情况。此外，当分为三个频带时，平均分类指标在大多数指标上表现出更好的性能。分类指标平均值的改善在分为三个频率频带时达到饱和。本文选择将其分为三个频率频带（[0.01, 0.04] Hz, [0.04, 0.07] Hz, [0.07, 0.1] Hz）来探索精神分裂症患者脑区信号异常与信号频率的关系。每个频率频带对应300个脑区特征，在拼接三个频带的特征后，每个受试者总共有900个脑区特征。同样，LASSO回归步骤重复1000次，每次回归总共保留了14个脑区特征：[0.01, 0.04] Hz频带10个，[0.04, 0.07] Hz频带2个，[0.07, 0.1] Hz频带2个。所有14个脑区特征在健康组和患者组之间均存在显著差异。
14个脑区中有11个出现在全频段筛选的异常脑区中，其余3个未出现的脑区特征（涉及脑区及编号如下：midf-97-'顶叶’，lowf-65-'后部’，higf-92-'颞极’）在健康组和患者组之间存在显著差异。即精神分裂症患者脑区65在[0.01, 0.04] Hz频带、脑区97在[0.04, 0.07] Hz频带、脑区92在[0.07, 0.1] Hz频带的平均信号强度存在异常。此外，为探索不同频段异常脑区信号的影响程度，还绘制了以网络为细分的14个脑区特征平均回归系数的条形图及带误差棒的折线图。属于控制网络和边缘网络的异常脑区在患者组中的平均振幅高于健康组；属于视觉网络和躯体运动网络的异常脑区在患者组中的平均振幅低于健康组。
3.4 使用DMD获得的模式重建和预测BOLD信号的有效性
用于执行DMD的BOLD信号值在[-50, 50]区间内，但集中在[-10, 10]区间。如图所示，使用每个受试者分解BOLD信号得到的所有139个模式进行重建的整体均方根误差（RMSE）平均为0.0044，而使用频率在[0.01, 0.1] Hz区间内的62/64个模式进行重建的整体RMSE平均为0.0068。良好的重建结果表明，在分解的时间窗口内，公式3可以明确表达每个时间点大脑BOLD信号活动的变化。
预测均方根误差序列在去除异常值后，用不同颜色的带误差棒的线进行可视化。当使用前130–139个时间点的时间序列进行预测时，每个剩余时间点的预测均方根误差保持在0.2以下；当使用前131–139个时间点的时间序列进行预测时，每个剩余时间点的预测均方根误差保持在0.1以下。使用100–139个时间点预测剩余时间点的总体效果是，前5个剩余时间点的预测非常有效，平均均方根误差低于0.3，从第6个时间点开始，预测效果越来越差，平均均方根误差呈指数增长。使用DMD提取的原始特征可以稳定预测接下来5个未知时间点（10秒）的BOLD信号活动变化，并且分解中使用的时间点越多，整体预测效果越好。
4 讨论
4.1 结果的可重复性及与现有结果的相关性
为检验结果的可重复性，还使用了‘Schaefer2018_100Parcels_7Networks_order_FSLMNI152_2mm’模板定义感兴趣区域，并获得了一致的结果。
精神分裂症以认知、情感和行为失调为特征，常伴有认知和社会功能障碍。功能研究发现患者脑区内激活异常和脑区间连接异常，但这些可能因患者表现出阴性或阳性症状而异。本研究中识别出的异常脑区可能在脑网络水平上代表了上述功能异常。与现有研究的进一步比较表明，平均振幅高于健康对照组的脑区可能对应于患者网络内功能连接增强或相应脑区激活水平增加。这为阐明精神分裂症脑功能改变的神经机制提供了新视角。
在本研究中，全频段的异常脑区在显著性/腹侧注意网络、边缘网络、视觉网络和躯体运动网络中具有相同的偏差方向；三个频段的异常脑区在控制网络、边缘网络、视觉网络和躯体运动网络中具有相同的偏差方向。这些空间一致的异常模式为进一步理解与精神分裂症相关的脑功能障碍提供了关键线索。
基于这些空间一致的偏差，我们可以辨别出精神分裂症患者大脑活动中一个更基本的系统级规律：精神分裂症患者在与内部动机、情绪和自我参照处理相关的网络（如显著性/腹侧注意网络、边缘网络和控制网络）中 consistently 表现出活动增强。相反，负责感知外部世界并对外部世界采取行动的网络（如视觉网络和躯体运动网络）则表现出活动减弱。这种“内高外低”的功能分离模式为理解该疾病的核心病理生理机制提供了关键线索。这不是随机的连接异常，而很可能反映了大脑整合内部状态与外部现实方面基本的协同功能障碍。这一发现与该领域的流行理论产生共鸣：首先，显著性网络（信息过滤的中心枢纽）的活动过度支持了“显著性网络功能障碍”假说，可能代表了驱动个体将内部想法误解为外部刺激，从而引发阳性症状的初始环节。其次，内部导向和外部导向网络集群之间的这种系统性失衡为精神分裂症的“脑网络失连接”理论提供了宏观层面的佐证。该理论认为，大规模脑网络之间功能整合和失连接机制的破坏是导致从思维障碍到精神运动障碍的广泛临床症状谱的基础。因此，本研究中识别的特征模式揭示了精神分裂症中一种潜在的跨网络系统性功能失连接模式，为其复杂症状谱的统一解释提供了新的神经生物学视角。
此外，患者脑区的异常信号具有频率依赖性，全频段的平均会丢失这部分信息。适当加入频率信息可以显著改善分类效果，并且还能发现患者隐藏在不同频段间的特征性异常信号脑区。
4.2 DMD与其他特征提取方法的比较
在生物医学信号分析中，自适应分解技术（包括基于小波的方法（如经验小波变换，EWT）和经验模式分解（EMD））已被证明能有效从各种信号中提取具有判别力的时频特征用于自动诊断。将这种数据驱动特征提取的范式扩展到功能神经影像，本研究采用的DMD具有分析大脑活动的独特优势：它从根本上设计用于表征时空动态。与上述技术生成的纯时频表示不同，DMD提取相干时空模式，直接捕捉大规模神经协调的演化模式，使其特别适用于研究如精神分裂症等疾病中改变的脑动力学。
在fMRI数据背景下，诸如PCA和ICA等方法已被用于识别形成功能连接的脑网络。然而，这些方法是静态的，因为它们将多元时间序列的连续时间点视为独立观测值。相反，DMD算法是动态的，结合了当今使用的两种最强大的数据分析工具的成熟优势：时间上的功率谱分析和空间上的主成分分析，这允许更清晰地区分在时间和空间上混合的模式，从而更好地提取动态数据的特征。动态分析发现的一些功能连接变化在静态连接分析中并未发现。此外，PCA和ICA都将信号视为特征时间序列的线性组合，这与大脑固有的非线性神经动力学不匹配。尽管识别DMD模式和特征值的过程也是纯线性的，但系统本身可以是非线性的，Koopman算子与DMD之间的联系可以证明非线性系统可以由一组模式和特征值来描述。因此，DMD提供了一个新的视角，使我们能够更深入地理解和分析脑活动的复杂动力学，这对于揭示脑功能的深层机制以及改善精神疾病的诊断和治疗具有重要意义。
为验证使用DMD提取的脑区平均振幅作为特征的有效性，采用了几种特征提取方法进行比较分析。这些附加特征包括基于DMD的特征的Gram矩阵（DMD-Gram）、基于PCA的特征的Gram矩阵（PCA-Gram）、基于ICA的特征的Gram矩阵（ICA-Gram）以及空间节点DMD特征（snDM）。需要注意的是，由于每个分量的地位平等，Gram矩阵不能与LASSO回归一起使用。在相同的训练/测试集划分下，使用线性SVM进行分类，并进行了三轮十折交叉验证。最后，进行了1000次置换检验以评估分类性能的显著性。结果表明，DMD-LASSO的分类性能与DMD-Gram、PCA-Gram或snDM-LASSO无显著差异，而所有四种方法（DMD-LASSO、DMD-Gram、PCA-Gram和snDM-LASSO）均显著优于ICA-Gram。
4.3 数据重建与预测
本文探讨了DMD算法在重建分解信号和预测未知信号方面的优良特性。DMD算法分解信号并给出模式、特征值以及每个时间点信号的显式表达式。输入不同时间点的值将给出该时间点的BOLD信号值，这使得能够重建分解后的信号并预测未知信号。原始BOLD信号数据包含一些噪声，使用特定频带的模式重建数据也可能是对BOLD信号去噪的有效方法。使用滤波后的模式可以更稳定地预测未来5个时间点的未知信号，并可能为时间点数据缺失的fMRI数据集提供更准确的插值。这进一步说明了大脑的复杂性和非线性，即DMD构建的动态模型仅捕获了分解窗口外几个时间步长内的脑活动。
4.4 DMD模式的序列相关矩阵与功能连接矩阵高度相似
在本研究中，进一步发现从每个受试者分解得到的所有模式之间的序列相关矩阵（phiC）近似等于其原始信号的功能连接（FC）矩阵。具体而言，每个受试者所有模式形成的矩阵维度为300 × 62/64，通过计算每个脑区中模式序列的皮尔逊相关系数，得到一个300 × 300的模式序列相关矩阵；而原始BOLD信号时间序列的相关矩阵则形成
