osl-dynamics：用于大脑快速动态活动建模的工具箱

1. 摘要

神经活动包含与认知相对应的时空结构。包括跨越大脑区域网络的振荡爆发和动态活动，所有这些都可以在数十毫秒的时间尺度上发生。虽然这些过程可以通过大脑记录和成像来获取，但由于其快速和短暂的性质，对它们进行建模在方法学上具有挑战。此外，认知事件的确切时间和持续时间通常是先验未知的。在这里，我们介绍 OHBA 软件库动力学工具箱 ( osl-dynamics )，这是一个基于 Python 的软件包，可以在快至数十毫秒的时间尺度上识别和描述功能神经影像数据中的循环动力学。其核心是机器学习生成模型，能够适应数据并在几乎没有假设的情况下学习大脑活动的时间以及空间和光谱特征。osl-dynamics融合了最先进的方法，这些方法可以并且已经被用来阐明各种数据类型的大脑动力学，包括磁/脑电图、功能性磁共振成像、侵入性局部场电位记录和皮层脑电图描记。它还提供了大脑动力学的新颖概括测量，可用于帮助我们理解认知、行为和疾病。我们希望osl-dynamics能够通过增强快速动态过程建模的能力，进一步加深我们对大脑功能的理解。

2. 引言

越来越多的证据表明神经振荡活动对大脑功能的重要性。神经振荡与认知过程有关，例如信息编码和处理以及注意力，并且在不同的意识状态下观察到不同的振荡活动。此外，神经振荡的同步已被提议作为一种沟通机制。神经振荡也是了解大脑功能障碍的有用工具。例如，在患病和健康人群的振荡活动中观察到变化。

神经振荡仍然有待充分理解的一个方面是其动态性质，特别是在快速时间尺度上。最近，有人提出，神经元群在 100 毫秒的时间尺度上表现出短时的振荡活动，而不是经典观点中持续振荡的幅度调制。这对于我们如何模拟认知和疾病中的振荡活动变化具有重要意义。不幸的是，可用于检测振荡数据中的方法有限，通常需要任意选择与频率内容、幅度和持续时间相关的参数。这些参数选择很大程度上影响了此类分析得出的结论。

此外，振荡爆发并不孤立于各个大脑区域。研究表明，诊断爆发可以发生在皮层网络中，并且网络中存在连贯活动的突发，在静息态和任务态中的持续时间约为 50-100 毫秒。对这些快速网络动态的精确了解是一种宝贵的见解，可以帮助我们理解认知过程；例如，特定功能性静息态网络的动态与记忆重放（记忆巩固中发生的≤ 50 毫秒的过程）有关。功能网络动态的变化也被证明可以预测行为特征和疾病。阻碍我们充分利用网络视角的关键障碍是动态功能网络的准确估计具有挑战性。这部分是由于有趣的认知事件的时间和持续时间以及功能网络中的相应活动尚不清楚。因此，我们需要依赖能够适应数据并自动识别网络何时激活的方法。

在这里，我们介绍 OHBA 软件库动力学工具箱 (osl-dynamics )，这是一个 Python 包，它解决了限制认知神经科学领域的两个影响深远的方法学挑战：突发检测和动态功能脑网络的识别。它通过部署数据驱动的生成模型来实现这一点，这些模型已被证明能够适应来自各种成像模式的数据，并且可以在很少的假设和快速的时间尺度上学习大脑活动的时空特征。

在突发检测应用中，osl-dynamics可以自动检测振荡突发，无需指定突发的频率、幅度阈值或持续时间。这使得osl-dynamics能够回答诸如以下的问题：振荡爆发何时发生；它们的振荡频率是多少；它们的特征是什么（例如平均寿命、间隔、发生和幅度）？在动态功能性大脑网络的检测中，osl-dynamics可以在快速时间尺度上自动检测网络动态，而无需做任何假设。这使得osl-dynamics能够回答以下问题：个人或群体表现出什么样的大规模功能网络；这些功能网络何时激活以及它们的特征动态是什么；哪些功能网络会激活以响应任务；个体的功能网络激活是否存在差异？除此之外，osl-dynamics可以使用与动态方法相同的方法，从更传统的静态（时间平均）角度来表征功能网络。

在这里，我们将使用公开的脑磁图 (MEG) 数据集来说明osl-dynamics的使用。然而，我们强调该工具箱的范围远远超出了 MEG，包含可以使用且已经使用的方法来阐明一系列应用中的网络和振荡动力学，包括脑电图、功能磁共振成像（fMRI），侵入性局部场电位记录和皮层脑电图。

3. 方法

3.1 生成模型

在动力学研究中，我们经常对时间序列的属性感兴趣，例如给定时间点的功率谱密度（PSD）、均值、协方差等。计算此值的常见启发式方法是使用滑动窗口。然而，这种方法仅利用感兴趣时间点周围的短窗口，并且在动力学的时间精度和属性的准确估计（通过足够大的窗口）之间进行权衡。过去的研究表明这种方法不足以研究功能连接的快速变化。在osl-dynamics中，我们采用基于生成模型的替代方法。这些是学习训练数据的概率分布的模型。在本报告中，我们将重点关注时间序列数据的两种生成模型：隐马尔可夫模型（HMM）和动态网络模式（DyNeMo）。这两种模型（下面进一步讨论）都在生成过程中纳入了潜在的动态潜在变量。训练期间的目标是学习最有可能生成观察到的时间序列的潜在变量（我们最小化变分自由能）。在此过程中，模型可以识别共享相同潜在变量的时间序列的非连续片段。汇集这些信息可以对数据的局部属性进行更可靠的估计。

HMM 的生成模型（如图 1A所示）的公式为：

其中θ t ∈ { 1, …, K}是时间 t 时的潜在状态，K是状态数，x t是生成的数据。p ( xt | θt )是观测模型。在这里，我们使用

其中µ k ∈ { µ 1 , …, µ K }是状态均值，D k ∈ { D 1 , …, D K }是状态协方差。时间序列中的动态是通过状态切换产生的，其特征是转移概率p ( θ t | θ t- 1 )。每个成对的状态转移形成转移概率矩阵，

osl-dynamics使用变分贝叶斯推理来学习最有可能生成观测数据的状态。这样做的优点是能够解释潜在状态的不确定性。有关osl-dynamics中 HMM 实现的更多信息，请参阅文档：https ://osl-dynamics.readthedocs.io/en/latest/models/hmm.html 。HMM已成功用于研究各种环境中神经影像数据的动态分析。

DyNeMo 是最近提出的模型，它克服了 HMM 的两个关键限制：互斥状态和有限的内存。DyNeMo 的生成模型（如图 1B所示）的公式为：

其中θ t是时间 t 的潜在向量（称为 logit ），x t是生成的数据。我们使用的观察模型是

其中µ j ∈ { µ 1 , …, µ J } 是众数平均值，D j ∈ { D 1 , …, D J }是众数协方差，J是众数，

是模式 j 的混合系数。潜在向量中的动力学是通过 p(θ t | θ 1: t- 1 ) 生成的，这是使用循环神经网络参数化的分布。具体来说，

其中 f 和 g 使用循环神经网络计算。osl-dynamics使用变分贝叶斯推理来学习最有可能生成观测数据的潜在向量。这是一种高效的推理方案，可扩展到大型数据集。有关osl-dynamics中 DyNeMo 实现的更多信息，请参阅文档：https ://osl-dynamics.readthedocs.io/en/latest/models/dynemo.html 。

经过训练后，两个模型都会揭示训练数据的动态潜在描述。对于 HMM，潜在描述是隐藏状态时间过程，它是在训练数据中每个时间点推断出的最可能的状态。对于 DyNeMo 来说，它是一个模态时间过程，是从训练数据推断出的每个模态的混合系数时间序列。我们将在第 2.5.1和2.5.2节中讨论如何使用这些潜在描述来总结训练数据中的动态。

图1：A) 隐马尔可夫模型 (HMM)。在这里，数据是使用隐藏状态 ( θ t ) 和观察模型生成的，在我们的例子中是由状态均值 ( µ k ) 和协方差 ( D k )参数化的多元正态分布。在给定时间点只能有一个状态处于活动状态。使用转移概率矩阵 (A ij ) 通过状态切换对动力学进行建模，该矩阵根据先前状态预测当前状态的概率。B）动态网络模式（DyNeMo）。这里，数据是使用模式（µ j和 D j）的线性组合生成的，并且使用循环神经网络（RNN：f 和 g）对动态进行建模，该网络根据特定混合比（α t）预测概率通过底层 logits ( θ t )计算。

3.2 数据集

我们使用两个公开可用的数据集：

CTF 静息态MEG数据集。其中包含使用 275 通道 CTF 扫描仪收集的静息状态（睁眼）MEG 数据。该数据集包含 65 名健康参与者的 5 分钟录音。它是在英国诺丁汉大学收集的，作为 MEGUK 合作伙伴关系的一部分。

Elekta 任务态 MEG 数据集。这包含在视觉感知任务期间记录的 MEG 数据。使用 Elekta Neuromag Vectorview 306 扫描仪记录了 19 名健康参与者的 6 个试次。该数据集是在英国剑桥大学收集的。

3.3 预处理和源码重建

图2显示了根据 MEG 记录估计源数据所涉及的步骤。这部分流程可以使用 OHBA 软件库(OSL)来执行，这是一个用于 M/EEG 分析的单独的 Python 包。应用于每个数据集的原始数据的确切步骤是：

1.MaxFilter（仅适用于 Elekta 数据集）。

2.带通滤波器 0.5-125 Hz。

3.陷波滤波器 50 Hz 和 100 Hz。

4.下采样至 250 Hz。

5.自动坏段删除和坏通道检测。

6.使用 EOG/ECG 通道的相关性来选择伪影组件，从而实现自动 ICA 清洁。

7.共同配准（使用 polhemus 头型点/基准点和结构 MRI）。

8.带通滤波器 1-45 Hz。

9.线性约束最小方差 (LCMV) 波束形成器。

10.分割到感兴趣的区域。在这项工作中，我们使用了 38 个区块。

11.对称正交化（纠正源泄漏）。

12.偶极子符号翻转（以对齐跨科目/运行的包裹时间课程的符号）。

13.下采样至 100 Hz（仅包含在突发检测管道中）。

这些预处理步骤适用于各种数据集。用于预处理和源重建的脚本可以在这里找到：https://github.com/OHBA-analysis/osl-dynamics/tree/main/examples/toolbox_paper。

图2. 首先，使用标准信号处理技术清理传感器级记录。这包括过滤、下采样和伪影去除。接下来，传感器级记录用于使用波束形成器来估计源活动。最后，我们对数据进行分割并进行校正（正交化和偶极子符号翻转）。这些步骤可以使用 OHBA 软件库执行：https: //github.com/OHBA-analysis/osl。

3.4 数据准备

我们通常在训练模型之前准备源数据。根据我们感兴趣研究数据的哪个方面，数据准备可能会有所不同。

振幅包络 (AE)如果我们对研究振荡振幅的动力学感兴趣，我们可以根据 AE 数据训练模型。在这里，我们通常对感兴趣的频率范围进行带通滤波，并使用希尔伯特变换的绝对值计算 AE。图 3B显示了当我们计算振荡数据的 AE 时会发生什么。我们可以看到 AE 数据跟踪振荡幅度的变化。

时延嵌入 (TDE)研究振荡的振幅动态并不能揭示不同区域如何通过相位同步相互作用的任何见解。为此，我们需要使用 TDE 准备数据。这通过包含原始通道的时间滞后版本的额外通道增强了时间序列。图 3C.I显示了一个示例。要执行 TDE，我们需要指定要添加的滞后通道数（嵌入数）以及要移动每个附加通道的滞后量。在osl-dynamics中，我们总是移动一个时间点，因此我们只需要指定滞后数。通过添加额外的通道，我们将原始数据的自相关函数（ACF）（以及互相关函数）嵌入到 TDE 数据的协方差矩阵中。图 3C.II对此进行了说明。我们在图 3C.III中绘制了取自 TDE 协方差矩阵的 ACF 和 PSD（使用傅立叶变换计算）。通过使用 TDE 数据，我们使协方差矩阵对原始数据中的振荡频率敏感。协方差矩阵还对通过非对角元素的跨通道相位同步敏感。对 TDE 数据的训练使我们能够研究通道之间振荡幅度和相位同步的动态。当我们准备 TDE 数据时，我们通常只对通过协方差矩阵寻找自相关函数/互相关函数的动态感兴趣，因此我们在生成模型中将均值固定为零。

有关在osl-dynamics中准备数据的更多详细信息和示例代码，请参阅教程：https://osl-dynamics.readthedocs.io/en/latest/tutorials_build/data_preparation.html。

图3. A.I）原始（模拟）时间序列数据。仅显示一小段（0.2 秒）。通道 1 (2) 是 15 Hz (30 Hz) 的调制正弦波，添加了N (0, 0.1) 噪声。A.II) 原始数据的协方差。B) 幅度包络 (AE) 数据（红色实线）和原始数据（蓝色虚线）。CI) 时滞嵌入式 (TDE) 时间序列。使用± 7滞后的嵌入窗口。C.II) TDE 数据的协方差。C.III) 使用 TDE 数据的协方差矩阵估计的原始数据的频谱特性。

3.5 一级和群体级分析

从源重建数据开始，我们通过两阶段过程研究数据集：

1)一级分析。在这里，我们的目标是估计特定主题的数量。在静息态（时间平均）分析中，我们直接从源数据计算这些数量。然而，如果我们正在进行动态分析，我们首先训练生成模型，例如 HMM 或 DyNeMo。请注意，HMM/DyNeMo 通常使用所有受试者的串联数据训练为组级模型。这意味着模型不知道数据来自不同的受试者，并允许模型汇集跨受试者的信息，这可以导致对动态量的更稳健的估计。我们使用模型提供的潜在描述和源数据来估计感兴趣的数量 - 这种方法称为双重估计。

2)集团层面的分析。为单个受试者估计的数量，例如网络指标或动态的汇总统计数据，用于对组进行建模。例如，这可以预测个体受试者的行为特征或特征，比较两组，或计算网络对任务响应的组平均值。通常，统计显着性测试是在组级别进行的，以验证任何观察到的差异或推断的关系不仅仅是偶然造成的。

我们将展示将五个管道应用于从2.2 节中提到的数据集计算的源重建数据的结果：基于 HMM 的突发检测管道（在2.5.1 节中讨论）；三个基于 HMM 和 DyNeMo 的动态网络分析管道（第 2.5.2 节中讨论）和一个静态网络分析管道（第 2.5.3 节中讨论）。HMM 和 DyNeMo 均已针对每个应用的模拟数据进行了验证。

3.5.1 突发检测

我们使用基于 HMM 的方法来检测振荡活动的爆发。在这种方法中，我们使用 TDE 准备源数据。典型的 TDE-HMM 突发检测流程如图 4所示。当根据训练数据推断 HMM 状态时间进程时，对特定状态的每次“访问”对应于突发或瞬态频谱事件，具有特定于该状态的频谱属性（例如β带功率的增加）。这种方法假设我们一次在单个通道（大脑区域）中寻找爆发；单独的 HMM 可用于检测每个通道中的突发。我们使用状态时间过程来计算表征突发动态的汇总统计数据。典型的汇总统计数据为：

平均寿命 ：这是状态处于活动状态的平均持续时间。

平均间隔：这是连续状态激活之间的平均持续时间。

爆数：这是平均一秒内某个状态激活的次数。

平均振幅：这是每个状态处于活动状态时源数据的 AE 的平均值。

图4. TDE-HMM 突发检测管道。这是在单个区域的宗地时间进程上运行的。为每个区域训练单独的 HMM。A) 通过执行时延嵌入和标准化（z 变换）来准备源重建数据。接下来，对 HMM 进行数据和统计数据的训练，这些数据和统计数据总结了根据推断的状态时间过程计算出的突发。B) 在组级分析中使用总结特定区域爆发的特定主题指标。

我们针对特定的州和主题计算其中的每一个。取所有状态激活的平均值。给定状态处于活动状态时，我们可以使用源数据来计算每种突发类型的 PSD。

3.5.2 识别功能网络

osl-dynamics为动态功能网络建模提供了更多选项。请注意，在这种情况下，我们对包含多个感兴趣区域（而不是单个区域）的活动的多变量数据进行训练，这就是我们在突发检测管道中所做的（第 2.5.1 节）。事实上，使用OSL 动力学对动态功能网络进行建模的一种观点是，它正在识别跨越多个大脑区域的突发。图 5显示了可用于动态网络分析管道的数据准备和生成模型的不同组合。我们将在下面讨论每个选项以及何时使用它们。

图5. 动态功能网络分析管道。A) 一级建模。这包括数据准备（如蓝色框中所示）、模型训练和事后分析（如红色框中所示）。第一级建模用于导出特定于主题的数量。B) 组级建模。这涉及使用第一级建模中的特定于主题的描述来对组进行建模。

AE-隐马尔可夫模型。如果我们有兴趣识别振幅的动态，我们可以使用 AE 数据进行训练。一旦我们训练了模型，我们就可以使用训练数据和推断的状态时间过程来估计主题和状态特定的网络（幅度图）。此外，我们可以计算汇总统计数据，以表征推断的状态时间过程的动态特征。这些汇总统计数据是：

占用率：这是每个状态处于活动状态的总时间的一部分。

平均寿命：这是状态处于活动状态的平均持续时间。

平均间隔：这是连续状态访问之间的平均持续时间。

切换率：这是特定状态每秒的激活次数。

我们针对特定的状态和主题计算其中的每一个。取所有状态激活的平均值。我们在第 3.2 节中展示了 Elekta 任务 MEG 数据集上的 AE-HMM 管道的结果。

时延估计-隐马尔可夫模型。我们可以使用 TDE 数据来研究相位同步的动态以及幅度的动态。在动态网络分析管道中，我们训练多元时间序列（即所有感兴趣区域的时间序列）。这意味着在 TDE 之后我们有非常大量的通道（嵌入数量乘以区域数量）。因此，我们经常需要进行主成分分析（PCA）来降维，以确保数据适合计算机内存。

在 TDE-HMM 管道中，我们可以计算与 AE-HMM 管道相同的汇总统计量。然而，为了估计功能网络，我们使用多维方法，我们使用源数据和推断的州时间进程来估计主题、地区和州特定的 PSD 以及交叉 PSD。然后使用 PSD 计算功率图并使用交叉 PSD 计算相干网络。请注意，我们还使用频谱分解方法来指定计算功率图和相干网络的频率范围。这涉及将非负矩阵分解应用于堆叠的主体和特定于状态的相干谱，以识别相干活动的公共频带。在本报告中，我们拟合了两个频谱分量，并且仅呈现第一个频段的网络，通常覆盖 1-25 Hz。我们将在第 3.3 节中看到在 CTF 休息和 Elekta 任务 MEG 数据集上应用 TDE-HMM 管道进行动态网络分析的结果。

TDE-DyNeMo。在此管道中，我们用 DyNeMo 替换 HMM 并在 TDE 数据上进行训练。与 HMM 提供的互斥状态描述不同，DyNeMo 推断模态时间过程，它描述了每个时间点每种模态的混合比例。这种混合描述使特定主题数量的计算变得复杂，例如网络和汇总统计。为了计算特定于模式和区域的 PSD，我们使用广义线性模型（GLM）的方法，其中我们将混合系数回归到（交叉）频谱图上。然后，我们使用模式 PSD 和交叉 PSD 分别计算功率图和相干网络。我们可以用平均值、标准差和成对皮尔逊相关性等量来总结每个模式时间过程的动态。或者，如果我们有兴趣计算与 HMM 相同的汇总统计数据（分数占用率、寿命、间隔、切换率），我们首先需要对模式时间过程进行二值化，这可以使用二分量高斯混合模型（GMM）来完成。请注意，与模式时间进程相关的另一个复杂问题是它不包含有关每个模式协方差的相对大小的任何信息。例如，如果模式协方差矩阵中的值相对较大，则混合比值较小的模式仍可能对瞬时协方差有较大贡献，我们通过重新规范化模态时间过程来解释这一点。我们在第 3.4 节中展示了 CTF 剩余 MEG 数据集上的 TDE-DyNeMo 管道的结果。

AE-DyNeMo。最后的选择是在 AE 数据上训练 DyNeMo。在这种情况下，使用 GLM 方法通过在滑动窗口 AE 时间过程上回归混合系数来计算幅度图。动态汇总统计数据的计算方式与 TDE-DyNeMo 管道相同。

当我们显示由上述每个管道推断的网络时，我们将对它们进行阈值设置，以仅显示最强的连接。在这项工作中，我们将使用数据驱动方法指定阈值，其中我们将二分量 GMM 拟合到每个网络中的连接分布。我们将其中一个分量解释为背景连接的分布，将另一个分量解释为异常的连接分布。

3.5.3 识别静态功能网络

osl-dynamics的一个特点是可以使用我们在动态方法中使用的相同方法来进行更传统的静态（时间平均）网络分析。这允许在静态和动态分析之间进行更直接的比较。为了对静态功能网络进行建模，我们只需要指定我们想要用来总结网络的指标，然后直接从源数据中计算这些指标。图 6显示了典型的静态网络分析流程。我们在第 3.5 节中介绍了 CTF 剩余 MEG 数据集上的静态网络分析流程的结果。请注意，对于静态网络，我们选择前 5% 的连接在每个图中显示，而不是我们用来对动态功能网络进行阈值处理的 GMM 方法。

图6. 静态功能网络分析管道。A) 源重建数据用于计算描述网络的度量。B) 特定于主题的度量用于对组进行建模。

3.6 运行间的差异

HMM 和 DyNeMo 通过最小化成本函数来训练（在osl-dynamics中，我们使用变分自由能）。通常，这种方法会遇到局部最优问题，即模型在训练期间可以收敛到数据的不同解释（潜在描述）。即，不同的状态/模式时间过程可以导致变分自由能的相似值。最终描述可能对更新模型参数和初始参数值的随机性敏感。

过去行之有效的处理此问题的策略是从头开始训练多个模型，并且仅分析变分自由能最低的模型。我们认为这个模型是对数据的最好描述，确保基于该模型的任何结论都可以在另一组独立运行的最佳模型中重现。

4. 示例分析

4.1 使用单区域 TDE-HMM 进行突发检测

图 4中的管道用于对左侧运动皮层中的单个包裹进行突发检测。源数据是使用 CTF 剩余 MEG 数据集计算的。所有受试者都在时间上串联并用于训练 TDE-HMM。结果如图7所示。我们从图 7A.I中的小波变换中看到，在这个时间序列中存在短暂的振荡活动。这说明使用传统的带通滤波和阈值方法来识别突发何时发生以及其中包含哪些频率是多么重要。我们使用三态 TDE-HMM 来以数据驱动的方式识别突发，而不是传统的突发检测方法。我们从推断的状态概率时间过程（图 7B.I）中看到，存在描述该数据的短暂状态。从图 7B.II中我们可以看到，每个状态都对应于独特的振荡活动。状态 1 被解释为非振荡背景状态，因为它的 PSD 中没有显示任何显着的峰值。状态 2 和 3 分别显示 5/θ频段 (1-7 Hz) 和 α/, β 频段 (7-30 Hz) 的振荡活动。图7B.III显示了每个状态概率时间过程与不同频段的AE的相关性（图7A.II）。基于此，我们将状态 2 识别为 5/θ突发状态，将状态 3 识别为 β 突发状态。从图 7B.IV中我们可以看到，这些突发具有从一百到几百毫秒不等的各种生命周期。

 

图7. 突发检测：单区域源重建的 MEG 数据（左运动皮层）显示短暂的振荡活动突发。A.I) 第一个受试者时间序列前 20 秒的动态光谱特性。A.II) 对 β 频带（上）、α 频带（中）和 5/ θ频带（下）上的时间序列进行带通滤波后计算出的幅度包络。B.I) 第一个受试者前 20 秒的推断状态概率时间过程。B.II) 每个状态的 PSD。B.III) 每个状态概率时间过程与不同频带的幅度包络的皮尔逊相关性。B.IV) 在受试者上的分布，用于表征突发的汇总统计。请注意，在计算平均幅度时，没有对源数据进行额外的带通滤波。

4.2 使用多区域AE-HMM检测网络动态

图 5中的 AE-HMM 管道应用于从 Elekta 任务 MEG 数据集中获取重建数据，以识别基于幅度的网络动态。所有受试者和运行都在时间上串联起来以训练模型。结果如图 8所示，其中包含对 HMM 状态时间过程进行组级分析的示例。

图8. 动态网络检测：在 Elekta 任务 MEG 数据集上训练的多区域 AE-HMM 揭示了与任务相关的具有快速动态的功能网络。A.I）对于每个状态，组平均幅度图相对于各状态的平均值（顶部）和绝对幅度包络相关网络（底部）。A.II) 说明第一科目前 8 秒的概率时间过程。A.III) 每个州的汇总统计数据的主题分布。BI）状态时间课程（维特比路径）围绕视觉刺激的呈现。水平条表示 p 值 < 0.05 的时间点。在排列测试中使用状态和时间点的最大统计池来控制族错误率。用于生成该图中结果的脚本位于：https: //github.com/OHBA-analysis/osl-dynamics/blob/main/examples/toolbox_paper/elekta_task/tde_hmm.py。

我们看到 AE-HMM 识别出具有快速动态的合理功能网络 ，通常寿命约为 50 毫秒（图 8A.III）。我们确定一个默认模式网络（状态1）；两个视觉网络（状态 2 和 6）；额颞叶网络（状态 3 和 7）；和感觉运动网络（状态 4）。

AE-HMM 以无监督的方式对连续源重建数据进行训练，即无需了解任务。HMM 训练后，我们可以围绕任务（呈现视觉刺激10）对推断的状态时间过程（Veterbi 路径）进行 epoch，并对试验进行平均。这给出了围绕视觉事件激活每个状态的概率。这如图8B.I所示。我们观察到在呈现预期的 50-100 毫秒后，视觉网络（状态 2 和 6）的激活显着增加（ p值 < 0.05）。我们还观察到视觉刺激后 300-900 毫秒额颞网络（状态 7）的显着激活（p值 < 0.05）以及视觉网络的失活（状态 2 和 6）。

4.3 使用多区域 TDE-HMM 检测网络动态

图 5中的 TDE-HMM 管道也应用于 Elekta 任务 MEG 数据集。所有受试者和运行都按时间连接并用于训练模型。结果如图9所示。

图9. 动态网络检测：在 Elekta 任务 MEG 数据集上训练的多区域 TDE-HMM 揭示了具有快速动态的光谱不同的功能网络。A.I）对于每个州，组平均功率图相对于各状态的平均值（上）、相干网络（中）和各区域的 PSD 平均值（下）、州特定的（彩色实线）和静态 PSD（即平均值）显示各州之间（黑色虚线）。A.II) 说明第一科目前 8 秒的概率时间过程。A.III) 每个州的汇总统计数据的主题分布。BI）状态时间课程（维特比路径）围绕视觉刺激的呈现。水平条表示 p 值 < 0.05 的时间点。在排列测试中使用状态和时间点的最大统计池来控制族错误率。用于生成该图中结果的脚本位于：https: //github.com/OHBA-analysis/osl-dynamics/blob/main/examples/toolbox_paper/elekta_task/tde_hmm.py。

对于高功率网络（状态 1-4），我们在 TDE-HMM 功率图（图 9A.I，顶部）和 AE-HMM 幅度图（图 8A.I）中看到相同的空间模式。我们可以从状态 PSD（图 9A.I，底部）看出，TDE-HMM 识别的网络表现出明显的光谱（振荡）活动。TDE-HMM 网络还具有快速动态（图 9A.III），寿命约为 50 ms。在图 9B.I中，我们可以看到我们能够重现使用 AE-HMM 所做的网络响应分析（图 8B.I）。

Elekta MEG 数据集是在视觉感知任务期间记录的。为了进行比较，我们对 CTF 剩余 MEG 数据集进行相同的分析。所有受试者都按时间连接并用于训练模型。图 10显示了将 TDE-HMM 管道应用于该数据集的结果。我们在休息时（图 10A）和任务中（图 9A）观察到类似的网络，这是功能磁共振成像研究的已知结果。我们还在图 9A.I和10A.I中包含了相干网络。我们观察到具有高功率激活的区域具有高连通性（一致性）。这些网络还具有快速动态（图 10A.III），寿命为 50-100 ms。

图10. 动态网络检测：在 CTF 剩余 MEG 数据集上训练的多区域 TDE-HMM 识别出与 Elekta 任务 MEG 数据集发现的功能网络相同的功能网络，并揭示了年轻组与老年组的动态差异。A.I）对于每个状态，组平均功率图相对于各状态的平均值（上）、绝对相干网络（中）和各区域的 PSD 平均值（下）、特定状态（彩色实线）和静态 PSD（即显示各州的平均值（黑色虚线）。A.II) 说明第一科目前 8 秒的概率时间过程并进行跑步。A.III) 每个州的汇总统计数据的主题分布。BI) 年轻组（18-34 岁）和老年组（34-60 岁）的汇总统计数据比较。星号表示 p 值 < 0.05。在排列测试中使用状态和指标的最大统计池来控制族错误率。用于生成该图中结果的脚本位于：https: //github.com/OHBA-analysis/osl-dynamics/blob/main/examples/toolbox_paper/ctf_rest/tde_hmm_networks.py。

为了说明我们可以从动态网络角度进行的群体层面分析，我们比较了两组：年轻组（18-34 岁）中的 27 名受试者和老年组（34-60 岁）中的 38 名受试者。图 10B.I显示了每组的汇总统计数据。我们看到老年组中感觉运动网络（状态 4）的占有率和切换率有所增加（p < 0.05）。老年组的视觉网络（状态 6）的平均寿命也有所缩短（p < 0.05）。对于默认模式网络（状态 1）和抑制状态（状态8）（p < 0.05），较旧的组的平均间隔分布也更广泛。

4.4 使用多区域TDE-DyNeMo的动态网络检测

图 5中的 TDE-DyNeMo 管道应用于 CTF 剩余 MEG 数据集。所有受试者都按时间连接并用于训练模型。结果如图11所示。请注意，对于 DyNeMo，我们发现学习 7 种模式（而不是 8 种）会带来更可重复的结果。因此，我们在图 11中展示了 7 种模式的拟合。

图11. 动态网络检测：在 CTF 剩余 MEG 数据集上训练的多区域 TDE-DyNeMo 揭示了光谱上不同的模式，这些模式比 HMM 状态更局部化并且在时间上重叠。A.I）对于每种模式，组平均功率图相对于各模式的平均值（上）、绝对相干网络（中）和各区域的 PSD 平均值（下）、模式特定的（彩色实线）和静态 PSD（即显示了各模式的平均值（黑色虚线）。A.II) 使用模态协方差的迹重新归一化模态时间过程（混合系数）。A.III) 通过连接每个受试者的时间序列计算出的重正化模式时间课程之间的皮尔逊相关性。A.IV) 重正化模式时间过程的汇总统计（平均值和标准差）的受试者分布。B.I) 年轻组（18-34 岁）和老年组（34-60 岁）的汇总统计数据比较。在排列测试中使用模式和度量的最大统计池来控制族错误率。用于生成该图中结果的脚本位于：https: //github.com/OHBA-analysis/osl-dynamics/blob/main/examples/toolbox_paper/ctf_rest/tde_dynemo_networks.py。

我们可以从功率图和相干网络（图 11A.I）中看到，与 TDE-HMM 相比，DyNeMo 识别出更多的局部功率激活和更清晰的网络结构。我们可以从 PSD（图 11A.I，底部）看出，这些网络还表现出独特的光谱特征。

从（重新归一化）模式时间过程（图 11A.II）中，我们看到 DyNeMo 提供的描述是混合比动态变化的重叠网络。这是对 HMM 提供的状态描述的补充。通过查看（重整化）模式时间过程之间的皮尔逊相关性可以理解每种模式的共同激活（图 11A.III）。我们观察到邻近区域的活动模式显示出更多的共同激活。我们使用图 11A.IV中的统计数据（平均值和标准差）总结了（重整化）模态时间过程。

为了在组级分析中比较 DyNeMo 与 HMM，我们使用 DyNeMo 特定的汇总统计数据（即重正化模式时间过程的平均值和标准差）重复年轻与年老的研究。图 11B.I显示了年轻（18-34 岁）和老年（34-60 岁）参与者的显着组间差异。我们可以看到感觉运动网络（模式 4）的模式贡献（平均重正化模式时间过程）增加，这反映了我们在 TDE-HMM 中看到的分数占用率的增加（图 10B.I）。我们看到，与p < 0.05 的 TDE-HMM 相比，DyNeMo 能够显示出更强的效应量，p < 0.01。DyNeMo 还显示时间网络（模式 5， p < 0.01 ）的变异性（重整化模式时间过程的标准偏差）有所减少。

4.5 估计静态功能网络

为了进行比较，我们还将典型的静态网络分析管道（包括静态功能连接）应用于 CTF 剩余 MEG 数据集。我们还考虑了年轻组与老年组级别分析中的静态视角与图 10中的 TDE-HMM 和图 11中的 TDE-DyNeMo提供的动态视角的比较，说明了能够进行静态和动态分析的好处在同一个工具箱内。

图 12显示了使用所有受试者计算的组平均 PSD (AI)、功率图 (A.II)、相干网络 (A.III) 和幅度包络相关 (AEC) 网络 (A.IV)。我们观察到 5 次幂在前部区域最强，α 次幂在后部区域最强。我们还观察到质量相似的相干性和 AEC 网络。特别是，我们在相干网络和 AEC 网络中看到 α 波段的枕骨连接性很强。

图12. 静态网络检测：osl-dynamics 还可用于执行静态网络分析（包括功能连接）。在 CTF 休息 MEG 数据集中，这揭示了年轻（18-34 岁）和老年（34-60 岁）参与者的静态功能网络中频率特定的差异。典型频带 ( δ、θ、α、β) 的组平均 PSD（对主题和地块进行平均，AI）功率图 (A.II) 相干网络 (A.III) 和 AEC 网络 (A.IV) 。B.I ) 老年人减去年轻人（上）和 p 值（下）的功效差异。仅显示至少有一个p < 0.05 的频段，其余频段均标有 ns（不显着）。B.II) 老年人减去年轻人的 AEC 差异仅显示 p  < 0.05 的边缘。在排列测试中使用频带和地块/边缘上的最大统计池来控制族错误率。用于生成该图中结果的脚本位于：https: //github.com/OHBA-analysis/osl-dynamics/blob/main/examples/toolbox_paper/ctf_rest/static_networks.py。

图 12B显示老年组（34-60 岁）减去年轻组（18-34 岁）的功率图 (B.I) 和 AEC 网络 (B.II) 之间存在显着差异（ p < 0.05）。我们观察到时间δ功率显着减少，感觉运动β功率增加。我们还观察到β带中感觉运动 AEC 显着增加（图 12B.II）。

我们在年轻参与者和老年参与者中看到的静态（时间平均）差异可以通过多种方式从静息状态网络的潜在动态中产生（图 10A.I和11A.I）。例如，静态功率的增加可能是由于特定网络的更频繁激活造成的。相反，网络的动态可能不受影响，并且网络内的功率可能会改变。使用 HMM 和/或 DyNeMo 研究动态网络视角有助于进一步了解静态差异是如何产生的。从 TDE-HMM 提供的动态网络角度来看，我们看到状态 4 的分数占用率有所增加（图 10B.I），这是一个在感觉运动区域具有高β功率和连接性（相干性）的网络。这与我们在此观察到的β功率和 AEC 连接性的静态增加一致；即随着年龄的增长，静态、 β-功率和连通性的增加可能与花在感觉运动网络上的时间的增加有关。TDE-DyNeMo 提供的视角显示，模式 4（代表感觉运动网络）的贡献随着年龄的增长而增加（图 11B.I ）。这是对静态、 β-功率和连接性增加的补充解释，因为感觉运动网络对老年参与者的整体大脑活动做出了更大的贡献。

5. 讨论

在第 3.1 节中，我们使用 TDE-HMM 以数据驱动的方式识别振荡突发，与基于幅度阈值的传统突发检测方法相比，假设要少得多。使用 TDE-HMM 等数据驱动方法的优点将在中进一步讨论。简而言之，使用传统方法，我们必须预先指定感兴趣的频率，并且我们可能会错过未达到任意阈值的振荡爆发。相比之下，TDE-HMM 对幅度不太敏感（能够更好地识别低幅度振荡突发），并且可以自动识别振荡频率。

在第 3.2 节和第 3.3节中，我们介绍了 HMM 在各种设置中识别的功能网络。这些网络是在快速（亚秒）时间尺度上从数据（无监督）中自动识别的，无需用户输入。我们发现了一组与任务（图 8）和人口统计（图 10 ）相关的合理网络。这些网络的可重复性非常好：跨多个 HMM 运行；跨越不同的数据准备技术（AE 和 TDE）；跨不同的实验范例（任务和休息）和不同的扫描仪（Elekta 和 CTF）。

鉴于我们使用 AE-HMM 和 TDE-HMM 观察到类似的网络（分别如图8和图 9），人们可能会问推荐使用哪种管道。TDE-HMM 方法能够对振荡幅度和相位同步的动态进行建模，而 AE-HMM 只能对幅度的动态进行建模。这意味着 TDE-HMM 通常是振荡动力学的更好模型。如果 TDE/PCA 数据训练的额外计算负载无法使用 TDE-HMM，则可能会首选 AE-HMM。

训练 HMM/DyNeMo 时必须做出的一个重要选择是状态/模式的数量。对于突发检测，我们通常感兴趣的是识别突发发生的时间点。这可以通过拟合两个状态 HMM 来实现：“开”和“关”状态。如果我们对多种突发类型感兴趣，我们可以增加状态数量。在这项工作中，我们选择了三状态 HMM 来保持接近开/关描述，同时允许多种突发类型。对于动态网络分析，我们需要少量的状态/模式来为我们提供数据的紧凑表示。常见的选择是在 6 到 12 之间。

osl-dynamics提供了两种用于检测网络动态的生成模型：HMM 和 DyNeMo。HMM 假设存在互斥的网络状态，而 DyNeMo 假设网络模式在每个时间点以不同的方式混合。虽然 DyNeMo 的假设可以说更现实，但 HMM 更强的假设具有简化分解的好处，这可以使解释网络动态更加直接。简而言之，HMM 和 DyNeMo 提供了网络动态的补充描述，其中任何一种都可能有用，具体取决于上下文。DyNeMo 确实具有通过使用深度循环网络来使用更丰富的时间正则化的额外优势。事实证明，这比 HMM能够捕获更长范围的时间结构，并且探索长范围时间结构的认知重要性是未来研究的一个有趣领域。可以使用下游任务来量化哪个模型更好。在这份手稿中，下游任务是诱发的网络反应和年轻与年老群体的差异。我们认为下游任务的更好性能表明模型更有用。

osl-dynamics还可用于计算静态网络描述，包括传统的静态功能连接。这使用与动态方法中的状态（或模式）特定网络估计相同的方法，使得动态和静态视角之间的比较更加直接。在3.5节中，我们使用这个特性将静态功能网络描述与动态视角联系起来。我们想强调的是，年轻人与老年人的研究被用作可以使用此工具箱进行的群体分析类型的示例，并且需要使用更大的人口数据集进行更严格的研究，以了解老龄化对功能的影响网络。第 3.5 节中的结果应被视为只是可能的老化影响的指示，可以在未来的研究中进行调查。在本报告中，我们重点介绍使此类研究成为可能所需的工具。

6. 结论

我们提出了一个用于研究时间序列数据的新工具箱：osl-dynamics。这是一个用 Python 编写的开源包。我们相信 Python 中这个包的可用性提高了这些工具的可访问性，特别是对于非技术用户而言。此外，它还不需要付费许可证。使用 Python 还使我们能够利用现代深度学习库（特别是TensorFlow ），这使我们能够将这些方法扩展到非常大的数据集，这是目前现有工具箱无法实现的。

osl-dynamics可以并且已经被用于广泛的应用和各种数据模式：电生理学、侵入性局部场电位、功能磁共振成像等。在这里，我们说明了它在突发应用中的使用使用 MEG 数据进行检测和动态网络分析。该软件包还允许用户在同一工具箱中研究时间序列的静态（时间平均）属性以及动态特性。OSL-动力学中包含的方法为动力学和群体层面的分析工具提供了新颖的汇总测量，可用于帮助我们理解认知、行为和疾病。

参考文献：osl-dynamics: A toolbox for modelling fast dynamic brain activity.