Nature Neuroscience：功能脑网络反映空间和时间的自相关性

神经科学中的高通量实验方法使测量复杂相互作用和多维模式的技术呈现爆炸式增长。然而，目前尚不清楚复杂的、反映涌现现象的指标能否归结到更简单的低维统计特征。为探讨这一问题，我们使用来自网络神经科学的复杂拓扑测量方法，分析了静息态功能磁共振成像（rs-fMRI）数据。在这里，我们表明空间和时间自相关是能够可靠解释众多网络拓扑指标的统计量。使用与被试匹配的空间和时间自相关构建的替代时间序列，几乎捕捉了这些拓扑指标中所有可靠的个体与区域差异。衰老过程中网络拓扑的变化由空间自相关驱动，而多种5-羟色胺类药物可因果地诱导相同的时间自相关的拓扑变化。这种对广泛使用的复杂性测量指标的还原式解释可能有助于将其与神经生物学联系起来。本文发表在Nature Neuroscience杂志。

正文

      随着神经科学数据变得更加复杂，用于分析这些数据的方法也变得更为复杂。复杂数据是否必须依赖复杂的方法？我们在这里关注的是功能连接（FC）的网络分析，即对静息态功能磁共振成像中大脑区域之间成对相关性的矩阵进行分析。FC是理解脑功能极为成功的工具。网络神经科学的非线性方法在功能连接的基础上，通过更高层次的抽象揭示了更多属性。这些方法中，功能连接通过图论被解释为网络，其中节点代表脑区，边代表这些脑区随时间活动的相关性。网络分析使得研究者能在更高层次的抽象水平上研究脑的拓扑结构，使用被称为“图指标”的统计量来研究如聚类、模块化组织以及脑全局活动的区域影响力等特征。

     网络神经科学深刻影响了我们对不同领域中rs-fMRI活动组织的理解。其中一个例子是衰老，这是一个复杂的过程，对脑的许多区域产生不同的影响。大量文献表明，网络拓扑反映了随着年龄增长在功能组织上的许多变化，包括全脑局部和整体连接的变化。除了网络拓扑的观察性和横断面差异之外，药理学的因果干预也显示出对被试内拓扑结构产生影响，这是转化应用的重要基础。例如，迷幻药物5-羟色胺受体激动剂麦角酸二乙酰胺（LSD）可以在因果意义上影响大脑网络的模块化组织和整体整合。然而，由于图论分析涉及多个抽象层次，尚不清楚个体内部或个体间的拓扑差异是否与更简单的时间序列特性差异有关。这可能导致通过基于MRI的图指标变化得出神经生物学结论时遇到困难。低层次因素，如静息态fMRI时间序列的基本统计特性，是否可能介导网络拓扑的变化？

     静息态功能磁共振成像（rs-fMRI）时间序列的基本统计特性在个体内和个体间存在差异。其中两个特性，空间自相关（SA）和时间自相关（TA），在神经影像数据中广泛存在，并且对数据统计分析具有显著影响。从定性上讲，SA体现了空间上邻近的数据点具有相似数值这一特征，而TA则体现了时间上邻近的数据点具有相似的数值。因此，SA和TA代表了变量之间的统计依赖关系，这种依赖关系的表现通常难以预测。尽管存在多种测量SA的方法，我们将其量化为空间中相似性随距离衰减的速率。同样，尽管存在多种测量TA的方法，我们采用时间序列中相邻时间点之间的相关系数作为测量方法。SA和TA之所以具有吸引力，是因为它们反映了众多的物理和生物学机制，包括脑的分子、结构、活动状态和组织属性。然而，尽管SA和TA在机制和组织层面具有重要性，以及其对统计数据分析的影响是已知的，目前尚不清楚SA和TA的个体差异是否能可靠地与更复杂的测量指标（例如网络神经科学的指标）的个体差异相对应。

     这里，我们表明两个简单的时间序列统计量——空间自相关（SA）和时间自相关（TA），能够解释网络拓扑结构中大部分的个体差异。我们证明了SA和TA的变化在不同受试者和不同脑区间是可靠的，并且SA与TA的个体差异与图指标的个体差异高度一致。

     为了理解SA、TA与图指标之间的关系，我们开发了一种仅以SA和TA为参数的时空替代时间序列模型。在将模型拟合于单个受试者的时间序列数据后，我们发现由受试者数据和模型生成的数据网络在图拓扑结构上具有高度相似性，证明了SA和TA对图论分析的直接影响。为了考察这种影响在实际数据分析中的意义，我们测量了随年龄变化而产生的网络拓扑变化。我们发现SA和TA不仅比传统图指标更敏感地捕捉了年龄相关的变化，而且我们的模型还能够将图指标的年龄相关变化分解为SA和TA的变化。此外，SA和TA对年龄相关的认知衰退也表现出敏感性。

     最后，我们展示了TA能够捕捉个体内部的变化，并且在超越图论框架的层面上与生物学过程存在因果关联。使用两种不同的5-羟色胺受体激动剂进行药理学的因果干预，结果显示出相同的TA下降的拓扑模式，而传统图指标并不能捕捉这一模式。此外，使用5-羟色胺受体拮抗剂则显示出相同拓扑模式下TA的增加，这表明相反的分子机制可导致相反的TA效应。因此，我们表明SA和TA是rs-fMRI数据的重要特征，与功能网络关系密切，为复杂脑网络拓扑提供了一种还原主义的理解方式。

结果
图论拓扑度量的低维性

     我们发现，从rs-fMRI构建的网络中，许多图拓扑指标在个体间高度相关。我们分析了多个由不同研究团队使用多种方法获得的神经影像学数据集，重点研究了人类连接组计划（HCP）数据集，并在耶鲁重测信度（Yale-TRT）数据集和剑桥老龄与神经科学中心（Cam-CAN）数据集中进行了验证。与以往工作一致，我们通过阈值化功能连接矩阵，仅保留前10%的强连接，并确保没有节点与网络断开，以构建脑网络。我们使用了若干图指标来量化网络拓扑结构：同配性、全局效率、传递性、模块化、平均聚类系数和平均局部效率。此外，我们还考虑了两个节点图指标：度数和中心性（方法部分详述）。

     在所有数据集中，大多数图指标在个体间彼此高度相关。然而，由二值化网络获得的无权图指标，也与功能连接（FC）的平均值（mean-FC）、方差（var-FC）和峰度（kurt-FC）高度相关。这一点令人意外，因为功能连接的二值化过程理论上破坏了有关mean-FC、var-FC和kurt-FC的所有显式信息。这些观察结果在所有数据集中均是一致的。因此，必须存在某些未被观察到的潜在因素，同时影响功能连接的统计特性和非加权图的拓扑结构。下面我们证明SA和TA就是这样的因素。

图1 功能性脑网络反映空间和时间自相关性

    a. 描述连接组分析流程与图指标的示意图。

    b. HCP数据集中所有受试者图指标之间的相关性，使用Bonferroni法进行家族错误率（FWER）校正后的双侧P值。

    c. SA-λ 和 SA-∞（上图）以及全局 TA-Δ1（下图）计算方法的示意图。计算 SA-λ 和 SA-∞ 时，首先按距离（D）将功能连接（FC）数值分组，然后将曲线 SA-∞ + (1−SA-∞) exp(−D/SA-λ) 拟合到分组数据上。计算全局 TA-Δ1 时，计算时间序列与向后移动一个时间点后的时间序列之间的皮尔逊相关系数。

    d. 使用组内相关系数（ICC）量化的图指标的重复测量信度。误差条表示95%的置信区间。每个柱状图后方有三个字符表示显著性：第一个字符表示相比SA-λ的显著性，第二个表示相较SA-∞，第三个表示与全局TA-Δ1比较，其中 "#" 代表P < 0.01，"+" 代表P < 0.05，"−"表示P > 0.5（基于双侧Bootstrap重采样法）。内嵌散点图展示了两个示例会话中跨受试者的相关性。

   e. 跨受试者的图指标与 SA-λ、SA-∞ 或全局 TA-Δ1 的相关性。“SA-λ + SA-∞” 和 “all” 分别表示使用两项或三项交叉验证线性模型预测图指标（随机选取50%的受试者用于训练），使用334个不相关受试者以避免亲缘关系混杂因素，显著性使用Bonferroni法进行家族错误率（FWER）校正。

   f. 跨所有受试者平均的区域 TA-Δ1 的脑区分布图。

   g. 各脑区区域 TA-Δ1 信度的分布。

   h. 各脑区区域 TA-Δ1 的信度分布图。

   i. 使用指纹识别分析正确识别受试者的平均比例。图中点表示从四次扫描会话产生的六种可能重复测试组合中的识别性能。

   j. 每位受试者的区域 TA-Δ1 与节点图指标在各个脑区间的相关性。

   k. 每个脑区的区域 TA-Δ1（取自 f）与其信度（取自 h）的关系。

    对于所有子图，除非特别注明，否则受试者数量 n = 883，* 表示 P < 0.05，** 表示 P < 0.01。箱型图显示中位数、第一和第三四分位数及范围，出于可视化清晰性考虑隐藏了离群值。FWER 指家族错误率（family-wise error rate）。

空间自相关（SA）

    空间自相关（SA）在神经科学中普遍存在，但定义通常并不明确，指的是相邻区域比距离较远区域更相似的现象。为了检验空间自相关（SA）的重测信度（即SA在同一个受试者不同的扫描会话中保持稳定的程度），我们开发了一种方法，在个体水平上将SA分解为两个组成部分：功能连接（FC）随物理距离衰减的速率（SA-λ）和两个远距离脑区间的平均相关性（SA-∞）。我们的方法首先根据空间距离将FC数据分组，然后为每个受试者的FC与距离曲线拟合出最佳的SA-λ和SA-∞（图1c），其中使用欧氏距离来同时涵盖大脑两个半球及皮层下区域（方程式1）（具体方法参见Methods）。重测信度使用组内相关系数（ICC）量化，它比较同一受试者多次观察的变异与所有受试者间的变异，其中1表示完全可靠，0表示随机水平。历史上对ICC数值的定性描述各有不同，通常0.5至0.7之间被称为“中等”到“良好”。

    我们发现SA在不同受试者之间可靠并且与图指标显著相关。与图指标相比，SA-λ和SA-∞均表现出较高的重测信度（图1d）（SA-λ的重测信度高于所有图指标，P < 0.01；SA-∞的重测信度高于未加权图指标，P < 0.05，经Bootstrap重采样方法检验）。这种效应不能用扫描过程中的头动（考虑头动后偏相关均>0.65，P < 10^−10，SA-λ和SA-∞均成立）或全脑信号回归（GSR）（扩展数据图1）来解释，并且在运动和脑区划分大小混杂因素最少的人类连接组计划（HCP）数据集中表现最强（扩展数据图2）。SA-λ和SA-∞均与加权和非加权图指标高度相关（图1e）。为了测试SA-λ和SA-∞的组合影响是否超过单独一个因素，我们构建了一个包含SA-λ、SA-∞和一个常数项的线性模型，并随机选择50%的受试者进行训练，以拟合每个图指标。该模型显著预测了留出数据集中的图指标（图1e），并且此结果在所有数据集中均成立（扩展数据图1）。因此，SA在个体间表现可靠并且具有拓扑学信息意义。

时间自相关（TA）
     TA 描述了 rs-fMRI 时间序列随时间的平滑度或记忆性，并已知在不同脑区之间存在异质性差异，且会同时影响功能连接（FC）和图拓扑结构。我们对每位受试者的每个分区脑区，采用该脑区时间序列在相邻时间点之间的皮尔逊相关系数来量化 TA，即滞后 1 的时间自相关（TA-Δ1）（图 1c）。TA-Δ1 是一种对 TA 的简单非参数度量。虽然从理论上说，TA-Δ1 只衡量时间序列中相隔一个时间点的相关性，但在实际中，它能够有效测量更长时间尺度上的相关性，包括对长记忆动力学的参数估计（扩展数据图 3）。由于预处理方法和重复时间（TR）都可能影响 TA-Δ1，因此只能在同一个数据集中对 TA-Δ1 进行比较。TA-Δ1 可以在个体脑区层面或整脑层面进行测量，我们分别称之为“区域 TA-Δ1”和“全局 TA-Δ1”。我们观察到枕叶和顶叶区域的区域 TA-Δ1 最高，而边缘系统区域的区域 TA-Δ1 最低（图 1f）。

    TA 在跨受试者层面上（无论是整脑还是区域层面）都表现出较高的可靠性。在整脑层面，我们将每个受试者所有脑区的区域 TA-Δ1 进行平均，得到该受试者的全局 TA-Δ1。结果发现，与各种图指标相比，全局 TA-Δ1 的重测信度非常高（图 1d）。在区域层面，我们也计算了各个脑区在受试者间的可靠性，结果显示，区域 TA-Δ1 的中位可靠性要高于其他节点图指标（图 1g）。此外，区域 TA-Δ1 的可靠性在全脑范围内也存在异质性差异（图 1h）。

    这种可靠性和异质性暗示，区域 TA-Δ1 可能被用于在群体水平上识别个体。为识别受试者，我们使用了“指纹识别”分析方法。对于给定的测量指标（如区域 TA-Δ1），我们将每位受试者在各扫描会话中的数值与所有受试者全部会话的数值逐一对比，找到该指标在脑区层面皮尔逊相关系数最高的那次会话。随后，我们统计有多少配对（会话）来自同一受试者。结果发现，在 1,765 个会话的庞大集合中，利用区域 TA-Δ1 进行指纹识别可以在超过 62.4% 的受试者中，成功匹配到唯一对应的那次扫描会话；相比之下，使用更传统的基于 FC 的指纹识别分析成功率为 62.6%（图 1i）。尽管这一数值低于文献报告，但我们数据集中可能的匹配数目级别更大，且与既往发现一致。我们还将其与利用若干节点图指标做指纹识别的结果进行了比较，后者无法达到同样水平的表现（图 1i）；在其他数据集中也获得了类似结果（扩展数据图 1）。此外，我们还检验了这种表现是否源于分区边界的不同或分区内部一致性（regional homogeneity）的可靠性差异，即将区域 TA-Δ1 与分区的区域同质性进行比较。结果显示区域同质性与区域 TA-Δ1 之间有中度相关，而其指纹识别能力在不同数据集间差异很大，且在 HCP 数据集中无法识别受试者（扩展数据图 4）。这表明，区域 TA-Δ1 的可靠性足以在群体中区分个体。

     在个体和区域层面上，TA-Δ1 都与图拓扑结构高度相关。在个体层面，我们发现全局 TA-Δ1 与各种图指标之间存在显著相关（图 1e）。为检验其与 SA 的共同影响，我们构建了一个包含全局 TA-Δ1、SA-λ、SA-∞ 及常数项的线性模型，在随机选取的 50% 受试者上进行训练以预测各图指标。结果表明，这一模型能显著预测所有的图指标（图 1e），并且在其他数据集中也几乎适用于所有图指标（扩展数据图 1）。在区域层面上，我们对每位受试者分别计算了区域 TA-Δ1 与多种图指标之间的 Spearman 相关系数。结果发现，不论是加权还是未加权的节点图指标，都与区域 TA-Δ1 呈高相关（图 1j）。尤其值得注意的是，节点的度（即该节点与其他区域建立连接的总数）与区域 TA-Δ1 的中位相关系数达到 0.89（图 1j）；且这一关系并非由分区大小所致（当考虑分区大小后，度数与区域 TA-Δ1 的中位偏相关系数仍为 0.83）。值得一提的是，这意味着在不直接研究网络拓扑结构的情况下也能够找到网络枢纽节点。

多种影响 TA 的因素
    跨脑区的 TA 异质性可能由多种潜在因素所决定。其中一个因素是脑区本身的内在时间尺度（intrinsic time scale），它反映了该脑区信号随时间衰减的速率。在 fMRI 和电生理学中，内在时间尺度一直被视为大脑动力学的关键因素之一。如果内在时间尺度是导致 TA 异质性的主要驱动因素，那么时间尺度更长的脑区应当表现出更高的 TA。然而，此前发布的内在时间尺度脑图与我们针对区域 TA-Δ1（图 1f）所得到的脑图并不吻合，这表明另一些因素也可能在决定 TA 方面发挥重要作用。

     除了内在时间尺度之外，脑区的噪声水平也会影响其 TA。当信噪比降低或噪声水平升高时，TA 的数值会下降。更一般地说，如果给任意时间序列添加白噪声，都会使其 TA 趋近于 0（见补充信息中的推导）。因此，由于大脑不同区域受到不同噪声源和噪声量的影响，这些差异也可能影响 TA。在一个极端情况下，如果噪声是唯一决定 TA 的因素，那么所有脑区都会有相同且缓慢的内在时间尺度，而每个脑区不同水平的噪声就会决定该脑区的 TA。与“噪声会塑造区域 TA”这一假设一致的是，我们发现 TA-Δ1 最低的脑区在 TA-Δ1 的可靠性上也最低（图 1k）。基于这一逻辑，TA 可能通过改变脑区对配对区域共享方差的贡献来影响功能连接（FC）。

    在接下来的内容中，我们利用这一原理，构建了一个时空模型（spatiotemporal model）来生成含有区域异质性噪声的替代时间序列。我们将其与“结合 SA 的内在时间尺度（intrinsic time scale with SA）”模型进行对比，后者将内在时间尺度视为驱动 TA 的主要因素。这两个模型在实际运用中最大的区别在于：在时空模型中，当 TA 减小时（由于噪声增大），SA 在局部产生的影响也会相应减小；而在“结合 SA 的内在时间尺度”模型中，SA 在构造时即与 TA 独立。我们还把时空模型与仅包含 SA 或仅包含 TA 的模型进行比较，以检验 SA 与 TA 的相互作用是否重要。最后，我们将时空模型与图论文献中一些常见的空模型（null models）进行比对。

用于替代时间序列的时空模型

      我们设计了一个时空替代时间序列模型，该模型在针对个体受试者拟合时，利用 SA 和 TA 来捕捉网络拓扑的个体差异。该模型作用于分区后的时间序列层面，这意味着可以在分析流程的多个层级进行比较（图 1a）。我们的模型可概括为少数几个步骤（图 2a）。首先，我们生成具有长记忆（1/f^α 频谱）的时间序列，使它们整体上具有较高的 TA；同时，根据大脑的几何结构对这些时间序列进行空间嵌入，通过增大相邻区域之间的相关性引入 SA，同时保留原有的频谱分布。这一步由我们所提出的相关谱采样（correlated spectral sampling）算法来实现（方法部分详述）。接着，我们向时间序列添加与各分区特异化方差对应的无相关白噪声，从而在不同区域异质性地降低 TA。这样得到的时间序列便可用与 rs-fMRI 相同的方式进行分析。因此，在我们的模型中，所有图拓扑结构的变化都必须由 SA 和/或 TA 的变化所致。

图 2：一个时空模型可捕捉大脑连接组的拓扑结构

    a. 示意图展示了我们用于生成具有 SA 和区域异质性 TA 之替代时间序列的时空建模框架。
    b. 针对一个示例受试者，将 FC 矩阵的特征值在对数坐标系下作图（黑色），与该时空模型（红色）进行对比。
    c. 示意图展示我们使用的模型拟合统计量——Linʼs concordance（Lin一致系数）在示例数据上的表现。当数值接近 1 时，表示模型既与数据呈相关，又在数值上非常接近（灰色）；数值接近 0 时，表示要么没有相关（粉色），要么存在较大平方误差（棕色）。
    d. 展示原始数据（黑色）与各模型（彩色）对应的 FC 矩阵（上方）及其图表示（下方）的示例。图以力导向布局（force-directed layout）进行可视化，将在拓扑上彼此接近的节点放置在较近位置。
    e. 对每种模型，将其与真实数据之间的 Linʼs concordance 值绘制出来。柱状图代表在四次扫描会话中的平均值，散点则表示每次扫描会话的 Linʼs concordance。为便于对比，黑色标记了来自同一受试者不同扫描会话之间的 Linʼs concordance，其中散点表示各会话配对。
    f. 在双对数坐标系下绘制的数据（黑色）与各模型的度数分布，用于对比。

   g. 各模型与真实数据在节点指标上的 Linʼs concordance 分布情况，按脑区进行统计。箱线图显示中位数、第一和第三四分位数以及范围，出于可视化考虑隐藏了离群值，n = 883。

    我们模型包含两个在单个受试者水平进行拟合的 SA 参数，分别对应无噪声的 SA-λ（即 SA-λ_gen）和无噪声的 SA-∞（即 SA-∞_gen）。我们还使用受试者实际观测到的区域 TA-Δ1 来确定需要加入的噪声方差。通过将 SA-λ_gen 和 SA-∞_gen 与 FC 矩阵的特征值分布进行拟合（该特征值分布本身不是图指标，但可被模型很好地捕捉，见图 2b），可以保证所有图拓扑的变化都必须来自于 SA 和/或 TA 的变化。

模型表现

     我们使用替代时间序列来评估模型对于加权和无权图指标的刻画能力。理论上，模型拟合的评估可基于多种标准，例如：模型对方差的捕捉能力（例如采用 R²），或者模型对个体差异的捕捉能力（例如采用皮尔逊相关系数）（参见补充表 1）。在此，我们使用 Lin 一致系数（Lin’s concordance）来评估模型的拟合度，该标准更加严格：只有当模型同时捕捉了总体方差和个体差异时数值才会接近 1；当二者任一较差时，数值即为 0 或负值。图 2c 展示了在模拟数据上，将 Lin 一致系数、相关系数和 R² 进行示意比较的结果。来自一个示例受试者的 FC 矩阵与网络图示表明，在我们测试的多种模型中，只有该时空模型能再现数据中观察到的复杂且不对称的模式（图 2d）。

     我们的模型很好地捕捉了图拓扑结构的重要特征。在加权和无权图指标方面，模型与数据之间的 Lin 一致系数均很高，意味着它不仅能够匹配个体差异，也能准确刻画图指标的具体数值（图 2e）。模型与数据的拟合度接近我们通过同一受试者不同扫描会话计算出的最大可靠变异范围（即使用同一受试者两次独立扫描会话在 Lin 一致系数上的结果，图 2e）。这种高 Lin 一致系数对应于模型与数据在图指标散点图上对角线附近的紧密分布（图 3）。和预期相符，我们还验证了模型能够匹配原始受试者的 SA 和 TA（扩展数据图 5）。此外，模型也能很好地再现节点度分布（图 2f），并且在个体水平上与多个节点图指标显著相关（图 2g）。

图 3：所有模型在图指标上的模型-数据相关性

     对于每种模型，将每位受试者的实际图指标与由时空模型预测的图指标进行绘制。散点图反映了图 2e 中所总结的关于 HCP 数据集中 Lin 一致系数（Lin’s concordance）的相关性。插图中给出了 Spearman 相关（r_s）和 Lin 一致系数（Lin）。“*”表示 Spearman 相关在双尾检验中 P < 0.05，“**”表示 P < 0.01。

    为了测试图指标可能与时间序列特性建立关联的其他方式，我们拟合了若干额外的模型，这些模型在扩展数据图 6 中有示意说明。
    首先，我们检验了模型中是否需要同时包含 SA 和 TA。为此，我们将时空模型中 TA 设为 0（即“仅 SA”模型），以及将 SA-λ 设为极小值且 SA-∞ 为 0（即“仅 TA”模型），并将这两种情形与我们原先的时空模型进行比较。
    其次，我们设计了一个模型来检验 TA 是否由脑区的内在时间尺度（intrinsic time scale）与 SA 共同决定（即“内在时间尺度 + SA”模型）。该模型为每个脑区估计了粉噪声 1/f^α 的指数，同时为全脑估计了 SA-λ 和 SA-∞，并利用相关谱采样（Methods）来生成替代序列。
   第三，我们考察了一种通过随机化相位来在所有尺度上重构 TA（“相位随机化”）的方法，是否能够与我们的模型表现相当。同样，基于我们在 mean-FC 和 var-FC 相关结果的发现，我们测试了一个能够匹配这两个统计量的模型（“Zalesky 匹配”），在先前工作中，该模型已显示出更接近脑网络的拓扑结构。此外，我们测试了一个能够匹配网络度分布的模型（“边洗牌”），这是图理论文献中广受欢迎的做法。最后，为确认我们的高拟合度是否源于在损失函数中使用了特征值，我们还创造了一种“特征值替代”（eigensurrogate）方法来构建替代的 FC 矩阵。该方法可保留 FC 矩阵的特征值分布，但随机化其特征向量，从而在完美复制 FC 矩阵线性维度的同时，用于检验效果是否仅由维度限制造成的零假设。

     这些替代模型都未能像时空模型那样充分捕捉网络拓扑（图 2e 和 3 以及扩展数据图 5）或度分布（图 2f）。将欧氏距离替换为测地距离（geodesic distance）对结果并无影响（扩展数据图 7）。在所有其他数据集中，只有在 Yale-TRT 数据上，特征值替代模型相较时空模型呈现了更强的解释力（扩展数据图 8 及补充表 1）。时空模型在对 Yale-TRT 和 Cam-CAN 数据集的整体拟合优于其他模型，强调了其普适性。

将自相关性与图拓扑结构联系起来
    我们的结果表明 SA 和 TA 能预测图拓扑结构，但由于对阈值化图进行构建和分析的过程是非线性的，想要在数学上严格建立这种关系并不容易。相反，我们试图通过考虑 SA 和 TA 对单条边的影响来为这一结论提供直观解释。就 SA 而言，其影响相对直接：SA 会增加彼此邻近的脑区之间的平均相关性，从而提高它们之间存在边的概率。但是，TA 会如何影响拓扑结构并不那么直观。

    一个统计学视角可以解释区域 TA-Δ1 与度（degree）之间的高相关性（图 1j），并说明高 TA 为何会产生网络枢纽。度由超过二值化阈值的相关边数目决定（图 1a）。换言之，即便节点的期望相关值（即节点层面的 mean-FC）不高，只要它的方差（nodal var-FC）足够大，也能使多条相关值越过阈值。因此，高方差比高平均值更容易增大跨阈值的概率。于是，任何可提升节点方差的过程同时也应提升该节点的度。理论上，两条时间上自相关的时间序列在两两相关上会呈现更高的方差补充信息 1.5），这一关系也反映在我们的数据里（图 1e,j）。因此，就单个节点而言，TA 使节点度增大的途径之一在于增大 var-FC。如果脑区 A 与脑区 B、C 的相关性均很高，那么 B 与 C 彼此相关的概率也更大。这意味着高 TA-Δ1 的节点在网络中更倾向于彼此相连，形成一个更聚类化（clustered）的网络拓扑。

    我们使用一个具有两个参数的“经济型聚类”模型（economical clustering model）来验证该推理。此模型根据距离和聚类拓扑，以一定概率直接将节点相连（补充信息 2）。在先前的研究中，该模型已被证实能够再现脑网络的若干拓扑特征，且其构造网络的流程与 rs-fMRI 流水线并不相同。结果发现，在“经济型聚类”模型中，SA 与 TA 的变化分别对应于更多短边和更高聚类边的倾向（扩展数据图 9）。换言之，在一个公认可再现 rs-fMRI 网络拓扑特征的图层面生成模型中，那些决定网络拓扑主要维度的参数与 SA、TA 这两个维度紧密关联。因此，我们可以从图拓扑结构的角度直接阐释 SA 与 TA。

健康衰老中的 SA
     我们已展示如何使用替代模型将不同受试者图指标的差异追溯到他们用于生成图指标的时间序列在 SA 和 TA 上的差异。由此，我们便提出相反的问题：如果图指标在群体层面上发生变化，那么是 SA 或 TA 导致了这些图指标的变化吗？由于在衰老研究中已有大量图论分析的积累，我们探讨了与年龄相关的图指标变化，并发现只有当 SA-∞ 参数发生扰动（而并非 SA-λ 或全局 TA-Δ1 参数），才能引起与衰老相似的图指标变化。

     我们分析了 Cam-CAN 数据集，该数据集包含了 800 多名年龄在 18 岁至 90 岁之间的受试者所采集的横断面 rs-fMRI 数据。尽管样本在年龄上差异巨大，但 SA、TA 及图指标三者之间的关系并未改变（扩展数据图 1），而我们的替代模型在这个数据集上依旧有效（扩展数据图 8）。由于该数据集中运动量与年龄高度相关（扩展数据图 2），我们在所有分析中都使用了对运动量进行控制的偏相关方法。多个加权和未加权的图指标都与年龄相关（图 4），其中全局效率（global efficiency）、var-FC 以及 kurt-FC 的相关性显著。年龄与 SA-λ 呈正相关，与 SA-∞ 呈负相关，而与全局 TA-Δ1 不相关（图 4a）。

图 4：在衰老过程中，SA 将功能性连接组拓扑结构与神经生物学联系起来

    a. 对图指标与年龄的偏相关进行了控制头动（motion）的分析。星号表示偏相关在双侧检验中的显著性（P < 0.05，n = 652）。误差条为 95% 置信区间。
    b. 可以对时空模型施加扰动，以测试每个 SA 参数是否介导衰老效应。上方：受试者数据预测，随着年龄增加，全局效率（global efficiency）和 kurt-FC 会下降，而 var-FC 会上升。根据 SA-∞ 和 SA-λ 随年龄的变化方式，分别对时空模型施加扰动。中图显示如果与年龄相关的变化仅由 SA-∞ 引起，则每个图指标随年龄变化的平均预测量；下图显示如果仅由 SA-λ 引起，则每个图指标随年龄变化的平均预测量。误差条表示该模型运行十次后的标准误。对于 SA-∞ 而言，图指标变化的方向（正向或负向）与实际年龄变化方向一致，而对于 SA-λ 则不一致。
    c. 显示了每个脑区的区域 TA-Δ1 与年龄之间的偏相关（控制头动）。
    d. 表示在全脑范围内，区域 TA-Δ1 与年龄之间偏相关的平均值（控制头动，n = 646 位受试者）。箱线图显示偏相关分布的中位数、第一和第三四分位数以及剔除离群值后的范围。表示在 Bonferroni FWER 校正后 P < 0.01。
    e. 度量指标与 ACE-R 评估之间的偏相关分析，控制头动与年龄。上图：对 SA 与全局 TA 的偏相关。SA-λ 与 SA-∞ 在 Bonferroni FWER 校正后显著，全局 TA-Δ1 不显著。下图：在所有加权或未加权的图指标中，只有 var-FC 在 Bonferroni FWER 校正后仍显著。误差条代表 95% 置信区间。 表示 P < 0.01，* 表示 P < 0.05，双侧 Wilcoxon 符号秩检验并采用 Bonferroni FWER 校正。FWER，家族错误率。

     我们的时空模型可用于确定网络拓扑变化在多大程度上由 SA-λ 与 SA-∞ 所主导。由于 SA-λ 与 SA-∞ 都与年龄呈显著偏相关，我们通过“逆向”运行该模型来理解这两个参数之中哪一个在衰老对于图指标的影响中起到中介作用。具体做法是，将这两个参数按衰老所预测的方向分别加以扰动。如果某个图指标的升降方向与实证中随年龄所观测到的方向相同，则说明此被扰动的参数在衰老效应中起到中介作用。由于 SA-∞ 会随年龄增大而下降，我们通过减少 SA-∞_gen 来模拟由 SA-∞ 介导的衰老过程，结果显示所有与年龄显著相关的图指标都随之发生变化，且该变化方向与实际的衰老方向相吻合（图 4b）。我们也尝试通过增加 SA-λ_gen 来模拟由 SA-λ 介导的衰老，但结果中图指标变化幅度很小，且方向与实际年龄变化方向相反（图 4b）。这表明，衰老对网络结构的影响主要由 SA-∞ 的变化而非 SA-λ 的变化所驱动。换言之，衰老对图指标的影响主要来自远距离脑区间的基线相关性水平，而非短距离随距离衰减的速率。

     由于全局 TA 并未随年龄发生显著变化，我们进一步探讨在区域层面上是否也是如此。我们计算了年龄与全脑范围内的区域 TA-Δ1 的偏相关（图 4c）。结果显示，在额叶子网络中，区域 TA-Δ1 随年龄上升而下降，但在小脑区域中则随年龄上升（图 4d），提示不同脑区的年龄变化模式存在差异。该差异可能反映了沿前—后轴方向的年龄相关结构改变。尽管效应量较小，这些结果依然说明，随年龄而出现的 SA 与 TA 变化，与图指标的年龄相关变化具有一致性。

亚临床痴呆标志
     虽然 SA 和 TA 与网络拓扑结构密切相关，但它们在解释临床症状（超越网络拓扑本身）方面是否同样有用？我们考察了 SA 与 TA 与痴呆早期症状之间的关系，关注在健康衰老效应之外，SA 或 TA 是否能预测认知衰退。为评估认知功能，我们使用了 Addenbrookes 认知测验修订版（ACE-R），这是一套用于亚临床人群痴呆筛查的认知测试量表。由于痴呆是参与 Cam-CAN 研究的排除条件，我们并不期待在 SA 或 TA 与痴呆标记之间能发现显著关系。我们针对所有受试者，计算了其 ACE-R 得分与 SA-λ、SA-∞ 以及全局 TA-Δ1 的偏相关系数，控制了年龄和头动。令人意外的是，结果显示 SA-λ 与 ACE-R 之间存在显著负偏相关（r = −0.184, P < 10^−5），表明更广泛的 SA 与更低的认知功能相关（图 4e）。SA-∞ 同 ACE-R 也存在较弱的相关（r = 0.118, P = 0.002），而全局 TA-Δ1 与 ACE-R 之间并无显著关系（r = −0.067, P = 0.09）（图 4e）。在所有图指标中，只有加权图指标 var-FC（r = −0.13, P = 0.001）显著相关（其余 P > 0.1）。这一结果突显了相对于图指标而言，SA-λ 和 SA-∞ 在敏感度上的优势。另一个值得注意的是，与主要由 SA-∞ 驱动的健康衰老相对比，此处的 SA-λ 效应却在此类认知衰退标记中起到了更突出作用。这些发现提示着临床症状可能与 SA 和 TA 存在某种关联。

药理学干预
     最后，我们检验了 SA 和 TA 是否能以因果的方式反映神经生物学过程。从理论上说，TA 和 SA 也可能仅由噪声、形态学或其他可重复的伪迹所驱动，而非真正的大脑动力学差异。通过被试内的药理学研究，我们可以对急性变化及其因果机制得出结论。因此，在两项人体药理学 fMRI 实验中，我们测量了 5-羟色胺受体激动剂 LSD和裸盖菇素（psilocybin）对神经回路调控后所导致的 SA 和 TA 的改变。在这两项实验中，受试者分别在不同日期以双盲方式服用药物或安慰剂，并在每次给药后两个时间点（早期与晚期）进行 rs-fMRI。即使样本量仅相当于 HCP 的 3% 以下，我们仍然观测到在药物作用下，SA、TA 与图指标间的基本关系依然保持（扩展数据图 10）。

     我们发现 LSD 和裸盖菇素都在全皮层范围内显著降低了 TA。LSD 与裸盖菇素在早期和晚期时间点均降低了皮层的 TA，在区域平均与被试平均的 TA 上表现相同（Wilcoxon 秩和检验，在所有条件下 P < 10^−40）（图 5a）。被试层面的药物对全局 TA-Δ1 的影响并非源于扫描过程中头动差异，而在此样本量下，并未能检测到药物对 SA 与图指标的影响（扩展数据图 10）。

图 5：血清素能药物的给药导致 TA 的改变。

    a. 皮层图显示在早期（左）和晚期（右）时间点，药物组与安慰剂组之间区域 TA-Δ1 的平均差值。紫色表示药物条件下的区域 TA-Δ1 低于安慰剂条件（LSD n = 24，裸盖菇素 n = 23）。
   b. 不同脑区在被试间的区域 TA-Δ1 均值随药物处理而改变的情况。箱线图显示区域 TA-Δ1 均值在被试间分布的中位数、第一和第三四分位数及排除离群值后的范围。* 表示 P < 0.05，** 表示 P < 0.01，双侧 Wilcoxon 符号秩检验。
  c. 测量每种药物与时间点的“大脑皮层图（药物减去安慰剂）”之间的相似度，使用 Spearman 相关系数。* 表示 P < 0.05，** 表示 P < 0.01；双侧置换检验同时控制 SA^[11^] 并使用 Bonferroni–Holm 校正以计算家族错误率（FWER）。
   d. 针对所有被试，汇总所有“药物减去安慰剂”对比所得到的第一奇异向量（SV）。右上：在每种药物条件下，第一奇异向量对个体差异所解释的方差均值。右下：在每种药物条件下的奇异向量得分均值。箱线图显示中位数、第一和第三四分位数以及排除离群值后的范围。* 表示 P < 0.05，** 表示 P < 0.01，双侧 Wilcoxon 符号秩检验，n = 646。FWER 指家族错误率。

    LSD 和裸盖菇素的致幻效应主要被归因于 5-羟色胺（5-HT）受体，尤其是 5-HT2A 受体的激动。对于 LSD 和裸盖菇素而言，如果它们共享相同的 5-羟色胺能机制，那么它们在整个皮层范围内对区域 TA-Δ1 的变化模式应当相似。实际上，我们发现无论是哪种药物、哪个时间点，皮层分布图之间都存在显著正相关（图 5c）。

    由于上述两种 5-HT2A 激动剂都在皮层范围内表现出减少区域 TA-Δ1 的特定拓扑模式，我们推断，一种 5-HT2A 受体拮抗剂应当在相同的脑区模式上产生相反效应，即增加区域 TA-Δ1。在 LSD 实验中，受试者在第三次访问时先被给予了选择性 5-HT2A 拮抗剂氯苯哒嗪酮（ketanserin），然后再服用 LSD（LSD+Ket）。我们发现，在晚期时间点，氯苯哒嗪酮对 LSD 所造成的区域 TA-Δ1 降低有显著阻断作用（图 5a），而在早期时间点，整体呈现出区域 TA-Δ1 的广泛上升（图 5a）。此外，在 LSD 条件下下降最明显的脑区，在给予氯苯哒嗪酮后则出现最显著的上升（图 5b）。这些与时间相关的 TA 变化与氯苯哒嗪酮的药代动力学一致，该药在初次给药后血药浓度会相对快速地下降。

    我们利用这些实验数据，构建了血清素能药物对区域 TA-Δ1 调控的整体皮层图。为整合所有实验条件的信息，我们对所有受试者的“药物对照”对比结果进行奇异值分解（SVD），提取了第一奇异向量。结果显示，在三个实验的早期与晚期阶段，该单一图（图 5d，左）能够解释个体间大约 50% 的方差（图 5d，右上）。在该图中，LSD 和裸盖菇素条件（包括 LSD+Ket 的晚期阶段）均呈负权重，而 LSD+Ket 的早期阶段呈正权重（图 5d，右下），与它们之间的相关结构相吻合（图 5c）。这些发现证明 TA 对血清素能药物的因果药理干预较为敏感，说明 TA 能够反映神经生物学上具有实际意义的差异。

讨论

    在本研究中，我们展示了以 SA-λ、SA-∞ 和 TA-Δ1 为参数的空间自相关（SA）和时间自相关（TA）是 rs-fMRI 时间序列的高度可靠特性，既与网络拓扑结构密切相关，也对其具有预测作用。当我们利用一个与被试匹配 SA 和 TA 的模型来生成替代时间序列时，这些时间序列所产生的网络能够匹配被试的图指标。我们还使用该替代时间序列模型来追踪图指标的年龄相关变化，发现尽管图指标本身尚不足以检测该变化，但 SA 和 TA 与亚临床痴呆症状呈现相关。此外，我们对多种血清素能药物进行因果药理学干预，结果表明 TA 会在皮层内呈现特异且稳定的调节模式，并具有较大的被试内效应量。我们预计，类似的 SA 和 TA 与网络拓扑关联的现象也会出现在其他来源于不同病症、疾病和药理学状态的数据集中。SA 和 TA 具备高可靠性、可解释性、较大效应量、与临床相关以及对神经生物学过程敏感等特点，使其成为 fMRI 潜在生物标志物的有力候选者。

     是什么因素驱动 SA 和 TA？我们发现，SA 和 TA 均可受到生物学因素（例如衰老）、药理学因素（如血清素能药物）以及方法学因素（如运动、分区大小和噪声）的影响。这些发现与先前在衰老过程中脑信号可变性方面的研究一致，也与其他非血清素类神经调质（如去甲肾上腺素、多巴胺、乙酰胆碱）的研究一致；同时也与方法学研究相契合，这些研究关注 SA 和 TA 的统计考量以及空模型形成复杂网络的能力。对于愈发引人关注的多样皮层梯度研究而言，SA 尤为重要。此外，SA 和 TA 还受到局部环路连接和内在时间尺度等性质的影响。

     我们还需更深入地理解生理性干扰过程对 SA 和 TA 的影响，包括 TA 与信号可变性的联系，尤其在衰老相关背景下。血氧水平依赖信号（BOLD）中的血流动力学响应函数不仅会对时间维度，还会对空间维度产生滤波效应，从而在衰老或血清素能药物的作用下导致 SA 和 TA 的改变。由于解剖结构在 SA 和 TA 中扮演了支架（scaffold）角色，我们所观测到的测试-重测可靠性和指纹识别表现也可能体现脑形态学在个体之间的差异。类似地，不同被试间功能网络组织的差异会带来分区边界的可靠错配，这有助于解释我们观察到的区域 TA-Δ1 的可靠性。我们的结果还显示，与上述结构和功能差异相关的分区内部区域同质性（regional homogeneity）与区域 TA-Δ1 呈正相关，但对指纹识别表现的解释能力因数据集而异。在 HCP 数据集中，它对指纹识别毫无解释力，而在 Yale-TRT 数据集中，它的指纹识别表现甚至优于区域 TA-Δ1。这说明在不同数据集中对 TA-Δ1 影响最深的潜在因素有所不同，也可能解释 HCP 与 Yale-TRT 两个数据集中模型表现的差异。TA 与方法学因素（例如重复时间 TR）之间也存在复杂关系：一方面，短 TR 使得采样时间点更紧密，对平滑的 1/f 功率谱而言，这会导致 TA 的上升；另一方面，短 TR 也可能造成信号更嘈杂，导致 TA 的下降。除此之外，预处理策略的影响也需系统化研究，以便在使用不同 TR 的研究中比较 TA。

    从总体上看，SA 和 TA 会约束神经信号的维度，而这些维度上的限制或许能够解释它们与网络拓扑之间的关联。以往研究已经证明，低维子空间能够捕捉网络拓扑结构，我们也通过经济型聚类模型展示了 SA 和 TA 如何与这些维度相契合。这些维度可能映射到网络拓扑的其他方面——例如，高度骨干（high-degree backbone）以及类似晶格的背景结构。图指标则为这些维度提供了复杂的非线性映射，使得我们能够在网络间比较拓扑性质。特征值替代（eigensurrogate）模型检验了一个更宽泛的问题：是否可以通过另一个具备同样线性维度限制的 FC 矩阵来复现图指标。尽管在多数数据集中它不如基于时间序列的模型表现得好，但也表现出相当亮眼的成果，甚至在 Yale-TRT 数据集中超过了我们的模型，这凸显了约束维度在塑造网络拓扑方面的重要作用。对特征值替代方法（包括对特征向量进行限制）的变体研究，或许能帮助我们进一步理解个体间脑拓扑结构的差异。

     本研究采用了分区分析的方法，这在一定程度上限制了对 SA 的测量精度。我们用于计算 SA 和 TA 的方法对于体素或顶点水平的分析可良好扩展，但若分区数目过多，则在计算上对拟合时空模型是不可行的。值得注意的是，SA-∞ 在概念上与 mean-FC 相似，而 mean-FC 已被证实与衰老存在关联；它也类似于全局信号（global signal），但即便在进行了全局信号回归（GSR）之后仍能保持可靠性，因此可能代表了全局信号在空间上的不均一性。尽管如此，即便是在分区水平上的分析，其灵敏度仍足以揭示图论分析无法捕捉的特征。例如，在我们药理学实验中，被试数量不足 25 人时，图指标并无显著变化，但 TA 的变化却揭示了功能组织方面高度显著的差异。

     我们的研究强调了在考察图指标这一复杂属性时，应与更为简单的属性——如 SA 与 TA——一起进行。日后在使用图指标的研究中，需要有明确的假设来指导如何在 SA 和 TA 的背景下慎重解读图指标。我们的研究也表明，仅基于 rs-fMRI 网络得到的图指标并不能直接视为区域间通信或信息处理的“标记”。我们的结果或许可拓展到更多超越图理论的功能连接分析情境，包括静息态网络和任务相关变化。我们的时空替代模型，以及先前已发表的在空间或时间上受到约束的空模型，也能更广泛地应用于神经影像学中的高维统计分析。本研究同样提示，SA 与 TA 在衰老及血清素能水平（serotonergic tone）等生物学过程中具有高度信息价值。历史上，SA 与 TA 常常被视为需要纠正的混杂因素，但我们的结果表明，它们不应仅仅被当作混杂因素，而是应被视作对连接组（connectome）而言至关重要且能提供信息的属性。