从静息态和任务态的大数据中发现动态大脑网络

大脑活动是对感官输入的反应和自发处理的动态组合。因此，无论一个人是否专注于外部强加的任务，这种大脑活动都会不断变化。之前，我们介绍了一种分析方法，允许我们使用隐马尔可夫模型 (HMM) 将任务态或静息态大脑活动建模为不同大脑网络的动态序列，从而克服滑动窗口方法带来的许多限制。在这里，我们提出了进一步研究的方法，使 HMM 能够处理非常大量的数据，从而在一系列不同的数据集中推断出可重现和可解释的动态大脑网络，包括任务态、静息态、MEG 和 fMRI，可能包括数千个被试。我们预计，人类连接组计划和英国Biobank等计划生成的大型公开可用数据集，再加上可以在这种规模上工作的计算方法，将为我们对健康和疾病中大脑功能的理解带来突破。本文发表在NeuroImage杂志。

1. 简介

在时间和空间尺度范围内理解大脑活动时间动力学的性质是神经科学的一个重要挑战。在研究任务数据时，目的是发现任务引发的神经基础和大脑机制，为此尽可能全面地将测量数据的时间过程与行为联系起来。也就是说，我们对任务引发的动态感兴趣。在这种情况下，通常会考虑同一任务的多次重复，以期表征和解释与某些基线条件相关的差异。据推测，大脑在不同的时间尺度上以在线方式适应任务，我们希望以成像模式允许的尽可能高的时间分辨率捕捉这些变化。在研究静息数据时，大脑没有进行预定义的任务，大脑仍将动态处理信息，使其活动适应当前对环境的感知，并结合其自身自发活动的产物。那么，在这种情况下，我们有兴趣表征自发动力学。无论哪种情况，考虑到大脑的复杂性和深度整合性，能够描述不同时间尺度的全脑网络活动的时间轨迹对于理解认知的最终基础至关重要（Buzsaki 和 Draguhn，2004 年；Bressler，2010 年）。

用于描述任务态和静息态中大脑网络动态的最常见分析技术是使用滑动窗口（Wendling et al.,2009 年；Allen et al.,2014 年）。滑动窗口方法（以及基于它的方法）存在许多限制（Hindriks et al.,2016 年），这些限制可能会使结论不够有说服力。特别是，他们需要预先指定感兴趣的神经过程发生的时间尺度，即窗口的时间宽度。这个选择是至关重要的，是两个相互冲突的标准之间的权衡：太长的窗口会错过快速动态，而太短的窗口将没有足够的数据来提供可靠的网络估计。滑动窗口方法的替代方法是隐马尔可夫模型 (HMM)，它可用于将大脑活动描述为离散大脑状态的动态序列，每个状态都具有独特的网络活动模式，包括功能连接和/或频谱内容（Baker et al.,2014 年；Vidaurre et al.,2016 年；Vidaurre et al.,Woolrich）。HMM 可以应用于任务数据，以提供对大脑动力学的丰富描述；例如，通过以完全无监督的方式（即不了解任务）估计 HMM，然后事后将 HMM 状态序列与任务时序相关联，以揭示任务相关的功能连接动态（Vidaurre et al.,2016 年）。这种策略可以揭示速度太快而无法通过滑动窗口分析看到的与任务相关的过程（Vidaurre et al.,2016 年）。HMM 还可用于静息态数据，以捕获在人群中持续重复出现的准静息态活动（Baker et al.,2014 年）。这允许分析某些动态属性如何因主题而异，例如状态之间的转换概率或状态占用的差异（即在每个状态中花费了多少时间）。图 1 展示了静息态和任务态中的 HMM。

在 HMM 的背景下，增加数据量有助于获得关于大脑活动动态特性的更丰富、更可靠的结论。例如，在任务中，进行更多的试验将使我们能够更好地理解与任务相关的大脑活动的时间及其逐个试验的变异性（部分原因是噪音，但也归因于有趣的认知过程比如学习）。幸运的是，人类连接组计划 (HCP)（Van Essen et al.,2013 年）和英国Biobank（Sudlow et al.,2015 年）等项目最近产生了数量空前的高质量数据，受试者数量达到数百或数千。然而，组级 HMM（在时间上连接所有主题的数据上运行）在如此庞大的数据集上进行训练在计算上是昂贵的。如果我们使用更复杂的 HMM 观察模型（即为每个 HMM 状态生成的数据模型），例如使用多变量自回归模型（MAR），该模型旨在捕获电生理学中特定于状态的多区域频谱信息，则会加剧这个问题数据（Vidaurre et al.,2016 年）。

在本文中，我们提出了标准 HMM 的替代方法，该方法使用随机变分推理方法，可通过大大降低其计算成本来应用于非常大的神经影像数据集。该算法普遍适用于不同数据模态所需的HMM框架的不同实例化。为了对其他研究人员有用，实现该算法的 Matlab 工具箱已公开发布。在目前的工作中，我们展示了该方法在来自 HCP 的 820 名受试者的静息态 fMRI 数据、来自Biobank 的 5248 名受试者的静息态 fMRI 数据以及来自 HCP 的 52 名受试者的 MEG 任务数据（最后一个数据集）的性能。总而言之，我们将这些示例与模拟数据一起使用，以证明拥有一种可以很好地扩展到大量数据的合适的计算方法可以显着丰富我们对动态大脑活动的描述。

2. 方法

2.1. HMM 和变分推理

HMM 是一系列模型，可以使用离散数量的状态来描述数据的时间序列，所有状态都具有相同的概率分布，但每个都有不同的分布参数。因此，这些状态对应于在时间序列的不同部分重复出现的独特大脑活动模式。对于每个时间点 t，状态变量指示每个状态在该时刻处于活动状态的概率。图 1a 举例说明了该模型；最上面是状态时间过程，即状态变量随时间的值；在底部，有一些状态的图形描述，在这种情况下，根据它们的平均激活水平和功能连接进行参数化。这个通用框架有不同的实例化，这取决于观察模型分布的选择。HMM 最常见的变体使用多元高斯观察模型，通常通过其均值 μk 和协方差 Σk 来表征每个状态 k 的分布（Baker et al.,2014 年）。重要的是要记住我们对网络的定义不同于独立成分分析 (ICA) 等提供的激活图。HMM 高斯状态不仅包含一个激活图（由 μk 编码），还包含一个协方差矩阵，可以解释为功能连接的网络矩阵（如图 1a 所示）。HMM 的另一个版本是 HMM-MAR，其中观察模型是自回归模型，因此状态由其频谱特征定义和驱动（Vidaurre et al.,2016 年）。在这种情况下，活动量和连通性都被确定为频率的函数。第三种可能性是 HMM-AR（具有自回归模型的 HMM），其中不对跨渠道交互进行建模。

无论选择哪种观察模型分布，HMM 通常都包含状态描述、状态时间过程（确定每个状态在时间序列中每个时间点处于活动状态的概率）和状态之间的转移概率（即从每个状态转换到另一个状态的概率）。因为在这里我们在所有连接的主体数据集上运行 HMM，所以状态和转移概率是在组级别定义的；然而，状态时间进程对于每个主题都是特定的，也就是说，状态可以在每个主题的不同时刻激活。由于模型每个部分的概率分布都依赖于所有其他部分，因此没有可用的封闭形式的解决方案。假设模型中的某些简化的流行推理范式是变分贝叶斯（Wainwright 和 Jordan，2008 年），它起源于统计物理学领域，并且在此之前，起源于 18 世纪由欧拉等数学家。变分推理方法在 HMM 概率分布中引入了某些因式分解，这样我们就可以迭代不同的参数组，使剩余参数保持不变，从而减少计算负担。目标是最小化自由能，该量包括真实分布和分解分布之间的 Kullback-Leibler 散度以及分解分布的熵。在 HMM 的背景下计算自由能的方程式可以在其他地方找到（Vidaurre et al.,2016）。

2.2.处理大数据集的随机变分推理

高斯过程中观测模型分布的估计意味着每个状态的 Q-by-Q 矩阵的反演，其中 Q 是通道数或时间序列（例如脑区）；对于 MAR （多变量自回归模型）模型，它需要对 PQ-by-PQ 矩阵进行求逆，其中 P 是 MAR 的阶数。在标准的变分推理方法中，任何一种情况都需要将整个数据集加载到内存中。状态时间过程的估计（对于任何 HMM 类型的模型，无论是使用变分推理还是最大似然法）基于所谓的 Baum-Welch（也称为前向-后向）递归（Baum 等al., 1970)，在计算了每个时间点每个状态的可能性后，需要对数据进行完整的顺序前向传递，然后进行完整的向后传递。虽然这个过程可以跨被试并行，但当时间序列很长和/或被试数量很大时，它仍然会很耗时。因此，HMM 的标准变分推理对于大型数据集可能具有挑战性，因为 (i) 估计观察模型所需的内存和 (ii) 估计状态时间过程所花费的计算时间。在本文中，我们利用应用于变分推理（Hoffman et al.,2013 年）的随机优化原理（Robbins 和 Monro，1951 年）并将其应用于 HMM，以开发一种非常有效的优化算法，解决计算困难，允许在短短几个小时内对许多主题进行 HMM 估计，并且内存要求适中。

图 1. HMM在 rest (a) 和 task (b) 上的工作流程图。

在这两种情况下，HMM 都会估计所有受试者或试验共有的几个大脑网络（或状态），以及每个受试者的特定状态时间进程，指示每个状态何时处于活动状态。在任务中，我们可以计算锁定到行为事件的状态均值激活，从而产生状态诱发反应，这对应于受试者处于每个状态的试验比例的时间进程。

标准变分推理保证自由能在每次迭代中减少，并最终收敛。随机变分推理反而在每次迭代时执行嘈杂且计算成本低的更新。尽管这些偶尔会导致小的自由能增量，但它们通常会改进模型。经过多次迭代，它可以保证收敛（Robbins 和 Monro，1951 年；Hoffman et al.,2013 年）。在这种情况下，“noisy”意味着我们将每次迭代都基于 M 个主题的随机子集：这些是我们需要保留在内存中的内容，以估计我们所说的临时状态观察模型，并且我们需要计算状态时间课程。（重要的是，为了获得一个临时状态观察模型，我们必须计算它的参数，就像实际使用了 N 个对象一样，以便估计的属性模仿标准变分步骤的属性）。这样，我们就有了观测模型的临时估计，它（由于高斯和 MAR（多变量自回归模型） 分布的可加性）可以与当前估计线性组合以形成新的估计。这种组合由一些标量 ρ 参数化，使得：

其中 φnew、φold 和 φinterim 分别表示观测模型的新的、先前的和中间的后验分布，并且 ρ 随着算法的进行而减小，因此在迭代 c 时我们有：α 和 β 是一些“延迟”和“忘记”参数。通过保持状态时间过程的充分统计，Σt Pr(st) Pr(st-1)，其中 st 表示时间点 t 的隐藏状态，并且 Pr (st) 表示每个状态在 t 处于活动状态的概率分布。受试者的抽样完全是随机的；在这里，我们建议（随机地）提升那些历史上采样次数较少的主题。我们通过以下等式来做到这一点其中 wi 是选择主题 i 的非标准化概率，ri 是主题 i 在先前迭代中被选择的次数（缩放使得 mini(ri)=0）和 τ ≤1 是控制我们阻止主题的程度的参数在当前迭代中经常被选中的那些。众所周知，HMM 优化可能会受到局部最小值的影响。出于这个原因，尽管随机推理可以通过更新中的噪音帮助避免平坦的局部最小值 (Hoffman et al., 2013)，但初始化起着至关重要的作用，因为它可以使优化过程远离差的区域参数空间。因此，我们需要一种在计算上在时间和内存使用方面都负担得起的初始化机制。我们在这里提出的初始化策略提供了一个相当不错的解决方案，而且计算成本不高。简而言之，它包括对主题子集单独运行标准 HMM 推理，并使用匹配算法将结果组合成单个解决方案。补充材料中提供了详细说明。该算法总结如下：1. 初始化观测模型。2. 重复：a.使用如（等式（3））中计算的概率 wi 随机选择 M 个对象。b.计算 M 个被试的状态时间进程。c.使用 M 个被试计算临时状态概率分布。d.使用 (Eq. (1)) 执行状态概率分布的近似更新。e.执行状态转移概率的精确更新。f.使用 (Eq. (2)) 更新 ρ。g.计算自由能。直到自由能收敛。

在本文中，我们使用了 α=5、β=0.7 和 τ =0.9。精确选择 α、β 和 τ 的影响有限。更重要的是 M 值的选择，我们根据数据集（见下文）进行了更改。非常低的 M 值会产生负面影响，因为更新变得过于嘈杂；非常大的值在计算上可能成本很高，而且通常不会带来太多好处。因此，选择的值是一种权衡。

2.3.模拟数据

我们首先使用合成信号来证明所提出的随机推理方法的有效性。我们使用 HMM 作为生成模型生成了两类信号。对于高斯观测模型，我们模拟了 6 个状态，每个状态有 10 个区域和一个随机生成的协方差矩阵。对于每个受试者，为了验证随机算法在存在一些受试者间变异性时产生与非随机算法相似的结果，通过添加第二个随机矩阵对协方差矩阵应用小扰动，从相同分布中采样。另一方面，对于 MAR （多变量自回归模型）观测模型，我们使用了 3 个状态，每个状态有两个通道/区域。这些状态是从对任务-MEG HCP 数据（如下所述）的实际 HMM-MAR 估计中获得的，我们观察到其中两个状态（红色和蓝色）正在捕获与任务相关的神经动力学，并且另一个状态（绿色）对应于基线状态。对于这两个观察模型，我们模拟了 200 个受试者的时间曲线，每个受试者有 500 个样本。

2.4. fMRI 数据

我们使用来自 HCP 的820名受试者的静息态 fMRI 数据，所有健康成年人（年龄 22-35 岁，453 名女性）在 3 T Siemens上扫描。对于每个受试者，可获得四次 15 分钟的 fMRI 时间序列数据，时间分辨率为 0.73 秒，空间分辨率为 2 毫米等距。在使用 FSL 工具去除人工噪声后（Jenkinson et al.,2012），我们使用组空间 ICA 获得覆盖皮质表面和皮质下区域的 50 个成分的分割；然后，我们使用这种分割将 fMRI 数据投影到 50 维时间序列中；这样的时间序列被标准化，以便对于每个扫描、受试者和 ICA 组件，数据均值为 0，标准差为 1。

我们还使用了来自UK Biobank的静息态 fMRI 数据，这是一项大规模前瞻性流行病学研究，旨在识别与许多疾病相关的风险因素和早期生物标志物（Sudlow et al.,2015 年；Milleret al.,2016 年；Alfaro-Almagro et al.,1000）。我们使用了包含 5847 名静息状态受试者（招募时年龄在 40-69 岁之间）的第一个公开版本。数据是在 3 T Skyra 中采集的，使用 32 通道接收头线圈，多波段加速度为 8 倍，持续 6 分钟和 10 秒，时间分辨率为 0.735 秒，空间分辨率为 2.4 毫米等距。预处理类似于 HCP 数据，包括运动校正、强度归一化、高通时间滤波和 EPI 反扭曲（Alfaro-Almagro et al.,1000）。然后，我们使用 ICA（Beckmann 和 Smith，2004），然后进行自动分类和去除人工噪声成分（Griffanti et al.,2014）。然后使用前 5248 个受试者的清理数据运行 Group-ICA，产生 100 个组，其中 55 个被手动分类为非人工噪声并用于 HMM 分析。执行对偶回归分析 (Beckmann et al., 2009) 以获得相关的特定节点时间序列，用于输入 HMM 分析。

图 2. 高斯和 MAR（多变量自回归模型） 情况下模拟数据的结果。

(a) 使用非随机方法和随机方法推断的状态时间过程之间的平均相关性（左），以及使用随机方法和基本事实推断的状态时间过程之间的平均相关性（右），对于具有两种不同的高斯观察模型不同批量大小的扰动大小和 MAR 观测模型。(b) 一个主题的基本事实和推断状态时间进程，使用具有较高的主体间可变性的高斯观察模型。 (c) Ground truth模型相对于标准和随机推理估计的平均激活（高斯分布的平均值）和功能连通性（协方差矩阵的非对角线元素）；每种颜色代表不同的状态，如果显示均值则每个点代表一个区域，如果显示功能连通性则代表一对区域。(d) 使用 MAR 观察模型，一个主题的基本事实和推断状态时间课程。(e) Ground truth模型（连续线）和估计模型（不连续线）的每个通道（左和右）和每个状态（蓝色、红色和绿色）的功率估计；上图对应标准推理，下图对应随机推理）。

2.5. MEG 数据

我们使用了来自 HCP 的 52 个 MEG 受试者，这些受试者同时具有静息状态和运动任务数据（Larson-Priora et al.,2013 年）。静息状态数据包括三个session（每个 6 分钟，睁眼），检查的任务数据有两个session，每个session 14 分钟。在进行 HMM 分析之前，时间序列也被归一化，因此对于所有受试者，它们的均值等于零，标准差等于1。

3. 结果

3.1模拟数据

目的是使用模拟数据通过比较使用非随机和随机算法的模型推理来验证算法的性能与标准HMM推理一致。我们针对 100 个不同的随机状态时间进程运行（标准）非随机推理和随机推理。由于输出中状态的排序是随机的，我们根据状态时间过程之间的相关性来匹配估计状态和实况模型。一旦状态匹配，我们使用匹配状态之间的平均相关性作为相似性度量来比较估计与基本事实，并将标准推理与随机推理进行比较。图 2a 显示了对于三种不同的随机批量大小，随机算法和非随机算法（左）之间的平均相关性以及均值随机估计和基本事实之间的相关性（右）。一般来说，即使在一个主题的批量大小的限制下，我们发现随机算法仍然推断出与非随机算法密切相关的状态时间过程。然而，正如预期的那样，对于较大的批量大小，两种算法之间的相似性更高（在批量大小等于受试者数量的限制下，两种算法完全等效）。

图 3. 对来自 820 名 HCP 受试者的静息态 fMRI 数据的随机 HMM 推断结果，使用高斯分布来描述每个状态。

(a) 三个示例状态的平均激活。

(b) 最大占用率直方图，一种旨在检查 HMM 是否能够表征数据动态的度量。

(c) 算法不同运行过程中状态时间过程的相关性，表明结果在运行过程中是稳健且一致的。

(d) 激活图的相关性和从数据集的单独半分割获得的估计之间的功能连通性，在 5 个随机分割中取平均值，以及从相关性从低到高排序的状态。

(e) 占用率（定义为受试者在每个状态花费的总时间）和每个状态的停留时间分布（即在每次状态访问中花费的时间），反映了数据时间动态的基本方面。

(f) 转移概率矩阵，反映每对状态之间转移的概率。

图 2b 和 d 分别显示了单个主题的推断状态时间过程的代表性示例，分别针对高斯和 MAR（多变量自回归模型） 情况（这些示例的准确度对应于 100 次运行的中值）。高斯情况下的真实平均活动（状态高斯分布的平均值）和功能连通性（状态协方差矩阵的非对角线元素）在图 2 中针对标准和随机估计（批量大小等于 50）绘制。图 2c，其中每种颜色代表不同的状态，每个点代表一个区域（如果显示平均值）或一对区域（如果显示功能连接）。两个通道的基本事实（实线）与标准推理和随机推理估计（不连续线）的状态谱信息如图 2e 所示。这些结果一起表明，对于各种配置参数，标准算法和随机算法之间的推断是一致的，两种算法都能够合理地反映基本事实。

3.2. HCP 静息 fMRI 数据

我们对来自 HCP 的 820 名静息状态 fMRI 受试者使用随机推断，以获得 12 种准静态大脑连接状态。我们使用标准工作站（配备四个 Intel Xeon CPU E5-2643 0 3.30 GHz 处理器）运行该算法 5 次，平均运行时间为 221 分钟（最小值和最大值分别为 208 分钟和 225 分钟）。我们使用了 M 1⁄4 30。下面的一些结果来自随机算法的选定运行。然而，不同的运行相对相似（见下文）。通过这些结果，我们的目标是说明：(i) HMM 可以提供的关于大脑动力学的信息类型，以及 (ii) 随机推理产生有用且重要结果的能力。请注意，在这种情况下无法比较标准推理和随机推理，因为给定数据集的大小，标准推理在计算上是不可信的。

图 3a 显示了十二个中的三个的平均激活图网络（状态），表示感觉状态、默认模式网络 (DMN) 状态和视觉状态（其余如图 SI-2 所示）。关于时间信息，最基本的输出是状态分数占用，定义为每个受试者在每个大脑状态中花费的时间比例。一个有效捕捉被试内时间动态的模型（与只发现被试间差异的模型相反）将被期望具有特定于被试的分数占用，这样单个状态就不会“支配”整个被试。换句话说，为了使 HMM 可用于描述大脑动力学，我们期望每个受试者的时间由不同的状态共享。反映满足这一最低要求的统计数据是最大入选率；也就是说，对于每个主题或扫描session，时间序列中花费时间最长的状态占用了多少时间。图 3b 显示了具有最大部分占用率的直方图。大多数受试者的最大占用率低于 0.4，表明大多数受试者需要多个状态才能得到最佳描述。

另一个重要的检查是随机推理的单独运行的一致性如何。这很重要，因为如前所述，已知标准 HMM 估计在某种程度上取决于初始化，并且我们的算法引入了一个额外的随机因素。为了研究这一点，我们匹配不同 HMM 随机估计的状态（使用 Munkres (1957) 的匈牙利算法，应用于状态时间过程的相关性），并收集重新排序的估计之间状态时间过程的相关性（即对于来自两次不同运行的每对匹配状态，我们得到这些状态时间进程之间的相关性）。图 3c 显示这些相关性很高（通常超过 0.75），反映了所提出的随机推理方法的稳健性。

进一步的稳健性测试，也将说明潜在科学结论的可靠性，是将数据集分成两半并分别运行 HMM。在这里，我们使用了 5 种不同的半分割，计算了每个分割的两个 410 主题 HMM 估计之间的激活图和功能连通性（协方差矩阵的非对角线元素）之间的相关性。图 3d 显示了每个状态的这些相关性（5 个分割的平均值），状态从相关性高到低排序。大多数状态在平均激活和功能连接的两个半分裂估计之间显示出高度相关性，除了一个状态，它与激活图的相关性相对较低。然而，在两种估计中，这种状态的平均激活非常接近于零，协方差矩阵反而会在状态处于活动状态时捕获数据的大部分不同的状态特定特征。此外，图 SI-4a 检查了半分裂之间转换概率的稳定性，这也非常稳健。

在评估了随机推断的基本有效性之后，我们更详细地分析了其中一次运行的时间特征。图 3e 显示了每个状态的部分入住率以及各状态停留时间的分布（指的是国事访问的持续时间）。从这个图中，我们可以看出一些差异，即一些状态的访问量比其他状态多。但是，停留时间方面的差异很小。最后，图 3f 显示了转移概率矩阵，表明某些（状态对）转移比其他转移更有可能（Allen et al.,2014）。总而言之，这些结果表明，当 HMM 与随机推理算法相结合时，可以在大型 fMRI 数据集中对大脑动力学进行可重复的建模。

3.3. Biobank 静息 fMRI 数据

接下来，我们对来自 Biobank 静息状态数据集的 5248 个 fMRI 受试者使用随机算法，再次执行 5 次算法运行，M = 250。平均运行时间为 144 分钟（最小和最大为141 分钟和 148 分钟）。在这种情况下，基于半分裂和 12 个状态的测试（参见 HCP 静息 fMRI 数据部分）显示 12 个状态中的 8 个在半部分是稳健的（参见图 SI-5）。HCP 数据集允许 12 个可靠状态，这可能是因为数据质量更高，每个受试者的扫描时间更长。然后基于结果的可靠性，我们在此展示具有 8 个状态的模型。其中三个代表感觉运动、DMN 和视觉网络，如图 4a 所示（其余如图 SI-3 所示）。其余图（图4b-f）与HCP结果类比，结论相似。HMM 可以捕获大脑活动的动态（图 4b），并且随机推理在算法的不同运行（图 4c）和数据集的随机半分割之间（图 4d 和 SI）是稳健和一致的；然而，这种一致性不如 HCP 数据集中明显，这可能是由于 HCP 数据的质量非常高。与 HCP 一样，各状态的部分入住率存在一些差异。然而，在这种情况下，停留时间的差异比 HCP 大，但仍然不大。

图 4. 来自 5248 个英国Biobank subjects的静息态 fMRI 数据的随机HMM推断结果，使用高斯分布来描述每个状态。

值得注意的是，尽管存在一些相似之处，但 HCP 和 Biobank 状态之间存在某些差异。这可能不仅是由于预处理和数据特征的差异，而且还因为数据被投射到不同的空间：虽然 Biobank 的状态图是体积图，但 HCP 图指的是皮质表面。实际上，ICA 分解强调了两个数据集中的不同区域。为了说明这种差异，图 SI-6 显示了“可表示性”图，其中，对于每个体素（或灰度坐标），我们计算所有 ICA 分量的绝对值之和。这表明每个 ICA 分解代表了多少区域。运行 HMM 的独立组件的这些明显差异可能解释了 HMM 结果的大部分差异。

3.4.任务 MEG 数据

最后，我们使用来自 HCP 的 52 个 MEG 受试者测试随机推理算法在任务中的性能。在这种情况下，我们使用 5 阶 MAR （多变量自回归模型）观测模型来描述状态，这样，如上所述，分割将基于数据的频谱信息。鉴于任务相对简单，并且因为我们只有两个数据通道（大脑区域），我们将 HMM 限制为仅推断 4 种状态。同样，我们使用 M=8 运行算法 5 次，平均耗时 166 分钟（最小值和最大值分别为 90 分钟和 229 分钟）。请注意，数据集的大小和模型参数的数量几乎允许使用标准推理方法，尽管计算成本很高。随机算法可以相对容易地处理这些数据和模型，而且，虽然我们在这里只关注两个区域以允许与以前的工作进行直接比较（Vidaurre et al.,2016），但它有可能被扩展，例如，在更多的大脑区域工作（例如跨越整个大脑）。

我们接下来表明，与之前使用标准推理（Vidaurre et al.,2016 年）所展示的类似，HMM 可以捕获任务调制的活动差异。图 5a（左上）显示了按钮按下事件前后的平均状态时间进程（具有标准误差），这可以解释为“状态诱发反应”，即每个状态的平均概率（跨试验）是多少每个时间点都活跃。两种状态，蓝色和红色，捕获大部分与任务相关的动态。通过应用统计测试定量证实了这一点，我们在每个状态和时间点测试了该时间点的状态分数占用是否与试验的其余部分（排列测试、显著性水平）显著不同）。 这些测试的结果描绘在图 5a 的顶部。可以观察到，随着时间的推移表现出差异的主要是红色和蓝色状态，这表明它们是由任务调制的。请注意，与 Vidaurre 等人的研究不同，此处的按钮按压有节奏且更频繁，可能会在运动中产生皮质活动夹带。由于前一个按钮按下的神经活动“泄漏”到下一个试验中，这种效果被 HMM 表示为状态时间过程的振荡行为。

图 SI-4 从 HCP (a) 和 UK Biobank (b) 数据集的五对半分割运行中获得的转移概率矩阵的相关性，我们将每个数据集随机分成两半，并分别在每一半上运行 HMM。每种颜色代表五个半分裂中的一个，每个点对应一对状态之间的转移概率。对于五个半分裂中的每一个，使用匹配算法在两次运行之间匹配状态。

图 5a 还根据功率（顶部中心）和相干性（右上）显示了状态的频谱特性。正如预期的那样，相干性的估计比功率的估计更嘈杂。有趣的是，代表两个半球 β 波段强烈激活的红色状态并没有像预期的那样在按下按钮的同时被抑制。同样，这可能是由于最后一次按下按钮引起的共振活动，这在新运动开始时带来了 beta 活动的运动后增量。可以结合时间和频谱信息来生成活动的 HMM 正则化时频表示。这如图 5a（底部）所示。

图 5. 来自 52 个 MEG 受试者的任务 MEG 数据的随机 HMM 推理结果，模拟按钮按下任务（1.2 秒试验间隔）期间两个运动皮层的活动。为了描述每个状态，我们使用了可以捕获数据频谱特性的多元自回归模型 (MAR) 分布。

(a) HMM 在模型推理过程中对任务信息视而不见，可以发现反映底层任务动态的状态；相应地，状态诱发反应（左上面板）区分概率增加的事件相关成分（例如蓝色状态，代表诱发反应）和其他显示减少的成分（例如红色状态，代表事件相关的去同步化）。当结合每个状态的时间（左上图）和频谱信息（右上图）时，该方法提供了数据的时频描述（下图）。值得注意的是，HMM 捕捉任务的节奏（即每 1.2 秒手部动作）。 (b) 状态开始概率，表示对于每个状态和时间点，该状态在该时间点变得活跃的试验比例；每条细线代表一个受试者，粗线是组平均值；顶部显示了主导状态是否比任何其他状态具有显着更大的发病概率的统计显著性。 (c) 说明未执行任务的数据的诱发反应。 (d) 算法不同运行的状态时间过程的相关性，证明运行非常一致。 (e) 从数据集的单独半分割获得的估计之间的功率谱密度的相关性，在 5 个随机分割中取平均值，以及从相关性从低到高排序的状态。 (f) 四个状态的部分占用率和停留时间。 (g) 转移概率矩阵。

作为部分占用率的另一种观点，我们还在图 5b 中查看了状态的开始（即状态变为活跃的瞬间）如何与任务相关。状态的开始可以解释为一个点过程，每个状态出现一个事件，指示该状态变为活动状态的时间。因此，对于每个时间点和状态，起始概率度量反映了该状态在该时间点精确激活的试验百分比。每条细线对应一个受试者，粗线代表受试者的平均值。为了可视化目的，曲线被平滑了。在图 5b 的顶部，每个灰点表示该时间点的主导状态比任何其他状态（排列测试，显着性水平 0.01）具有显著更大的起始概率，表明它是红色和蓝色状态主要达到意义。该图给出了状态动力学的补充视角。例如，它显示蓝色状态平均而言在按下按钮的那一刻变得活跃，在 ~+0.25 秒时达到部分占用率的峰值（图 5a），并且稍晚于红色状态±0.5 秒。

图 5c 显示了从 MEG 静息状态数据中获得的 HMM 估计的平均状态分数占用（具有标准误差），其中没有执行任何移动，并且试验被人为分配给静息状态session。在这里，状态没有显示任何振荡模式或任何其他一致的时间变化（即平均状态时间进程在整个试验中是平坦的，并且在全局时间信息方面状态之间的唯一差异是整个平均分数占用率所有时间点）。最后，图 5d 和 e 表明，在算法的不同运行和数据集的随机半分割之间的估计非常一致（对于半分割，相关性是根据从 MAR（多变量自回归模型） 系数获得的功率谱密度计算的）。接下来显示其中一个运行的附加时间信息，显示分数占用和停留时间分布（图 5f）和转移概率矩阵（图 5g）。

补充：

图SI-6 可表示性图反映了 ICA 分解代表不同大脑区域的程度。对于每个灰度坐标（左）或体素（右），我们计算了所有 ICA 组件的绝对值之和，因此较高的值（较暖的颜色）表示该区域更好地表示。

4. 讨论

目前正在努力公开发布包含数百甚至数千个主题的大型数据集。此外，协议的改进和标准化，以及广泛数据共享的增加，可能使合并和融合来自不同记录方式的神经影像数据成为可能。更多数据的可用性允许进行更复杂的分析，这将使我们能够首次提出特定问题。由于这些原因，拥有可以扩展到大型数据集的方法至关重要。在本文中，我们提出了一种可以利用大数据集可用性的方法解决一个难题：任务态和静息态中动态连接的表征。

HMM 方法在过去已被证明是有用的（Baker et al.,2014 年；Vidaurre et al.,2016 年；Vidaurre et al.,Woolrich）。在静息状态下，全脑静息状态网络的动态特征已在 MEG（Baker et al.,2014 年）和 fMRI（Vidaurre et al.,Woolrich）中进行了表征，例如，揭示了状态或网络之间的转换很远来自随机。在任务中，HMM 能够找到具有特定频谱特征的短暂、短暂的状态，这些特征在简单的运动任务中有意义地表征了初级运动皮层动力学。在这里，我们大大增强了方法，以允许在任务态和静息态中从非常大到非常大的数据集表征动态大脑网络。

HMM 提供了对脑网络动力学的丰富描述，不受滑动窗口方法的限制。更具体地说，随着推断的 HMM 状态数量的增加，我们有可能捕获数据中包含的大部分信息，而不受预定义的时间尺度（窗口的宽度）的约束，并且没有任何估计问题滑动窗口方法（即仅基于几个数据点的估计的不确定性）。例如，随着 Biobank 数据集向 100,000 名受试者的目标增长，使用提议的随机推理方法，我们有信心估计具有更多状态的 HMM，同时仍然有足够多的数据来足够好地估计每个状态。将进一步提升数据潜力的两个有趣的发展可能是：（i）使用一些原则性定量方法组合来自不同数据集的模型，以及（ii）在单个中使用不同数据集估计联合模型的可能性推理过程。这将需要足够统一的预处理管道并将数据映射到一个共同的大脑空间（目前 HCP 和 UK Biobank 不是这种情况）。虽然我们的方法可以很容易地应用于静息态和任务态数据，但静息状态下功能动力学的研究并非没有争议。例如，Laumann 及其同事 (Laumann et al., 2016) 有理由认为，大多数观察到的功能连接性变化是由于认知内容的真正转变以外的因素造成的，即头部运动、睡眠和采样变异性。

总之，这里介绍的方法为基于 HMM 的新方法研究打开了大门。例如，许多不同模型的组合以获得比其部分表现更好的整体，这在监督学习中是一个成功的故事（Hastie et al.,2003 年）。在这里，在无监督学习场景中，我们可以利用算法的计算效率来获得大量的 HMM 集，这些 HMM 可以以某种方式组合以提供更丰富的数据描述。如何组合这些模型，以及我们可以通过使用单个 HMM（可能具有不同的参数化）的这种合并来研究数据的哪些新方面，将是未来研究的主题。