事件标记窗口通信：从神经时间序列推断活动传播

摘要

追踪活动或信号扰动如何在神经系统中传播，对于深入理解大脑区域间的通信机制至关重要。然而，现有的分析方法尚无法系统地从神经时间序列记录中推断区域间的活动传播。为此，本文提出了一种新的框架，即“事件标记窗口通信”(EWC)，通过追踪自发性内源区域扰动的统计结果，来推断神经元素之间的活动传播。EWC通过对神经时间序列进行子采样，并利用成熟的功能连接性度量方法来量化统计依赖关系，从而追踪这些扰动的下游效应。本研究通过神经动力学模拟测试了EWC，并展示了其在多种模型配置下都能准确还原方向性信号传递的真实模式。通过与传递熵等更先进的功能连接估计方法进行对比，本研究还证明了EWC在捕捉活动传播方面具有较高的计算效率。最后，本研究展示了EWC在脑磁图(MEG)记录中推断全脑活动传播图的效用。此外，使用EWC计算的网络与通过传递熵推断的网络具有高度相关性(r=0.81)。总的来说，EWC框架具有较高的灵活性，可应用于多种功能神经影像技术获取的活动时间序列，为神经通信研究开辟了新的途径。

引言

神经元素(如神经元、神经元群和灰质区域)之间的通信是大脑功能的核心。理解信号在网络化神经系统中的动态和灵活传输机制，仍然是当前研究面临的重要挑战和热点问题。该领域已提出大量脑通信模型，涵盖了从神经振荡的信息传递到基于连接组通信的网络度量等多种方式。尽管存在众多模型，但对这些模型的验证工作尚显不足，而且也不清楚哪些模型能够真实反映神经动力学和信号传递的模式。模型验证的一个主要障碍是目前无法从神经活动记录中推断信号传递事件。由于现有神经影像技术的限制，活动传播追踪通常仅用于微观神经元网络中的小规模受控刺激实验。虽然这些技术在全脑尺度上的应用仍然有限，但整合现有数据来研究大脑活动可能会为我们提供更深入和可靠的见解。

确定信号、区域活动或内源性/外源性扰动通过解剖基质在区域间传递的程度，可为区域间通信原理提供重要见解，并有助于验证现有的计算模型。简而言之，一个区域“X”如何与区域“Y”进行通信，理论上可以通过测量Y在X发生扰动后的活动来确定。将这一概念扩展到时间维度，X的一系列扰动可能引发Y的一系列耦合响应，从而形成区域间的统计依赖关系。功能连接性(FC)可作为评估这种依赖性的工具。从广义上讲，FC通过Pearson相关系数或互信息等度量来量化神经活动时间序列的统计依赖性，捕捉两个脑区在一定时间内的动态同步程度。早期FC研究主要在整个时间序列上进行估计，以揭示区域间显著且稳健的统计关联。该方法后来扩展到通过滑动窗口定义的扫描子样本估计FC，从而捕捉更小时空尺度上的统计关联，称为动态FC(dFC)，以区别于“静态”FC。动态和静态FC估计已广泛应用于模型构建、验证以及实验和临床研究。然而，在存在自发性、离散性和方向性信号传递事件的系统中，当采用静态或动态FC估计方法时，反映通信的统计依赖性可能会被与通信无关的信号段掩盖。信息论方法(如传递熵等)有望解决此问题，但在进行全脑网络推断时，通常需要大量的数据和计算资源。

本研究旨在开发一种高效的方法，通过将无向FC估计的实用性与活动传播追踪所提供的动态分辨率相结合，从而推断活动传播模式。本研究将该方法称为“事件标记窗口通信”(EWC)。具体而言，EWC包括识别神经记录中的显著事件(“刺激”或扰动)，随后在短时间有序窗口/子样本中估计信号的统计依赖性(FC)，有效地捕捉下游目标对这些事件的响应。通过这种方式，EWC捕获了活动传播的间接统计模拟，而非时空局部化的刺激-响应测量。本研究首先在嵌入真实信号传递模式的简单计算机模拟网络上验证EWC的有效性，并证明了该方法允许使用无向FC度量来捕捉方向性交互作用。然后，通过比较网络推断运行时间与节点数量的函数关系，研究该方法的计算可行性。最后，本研究展示了该方法在源定位静息态脑磁图(MEG)记录中的实际应用。

方法

数据集

通过ConnectomeDB平台，从人类连接组项目(HCP)中获得了30名受试者(年龄22-35岁，17名女性)的静息态脑磁图(MEG)扫描数据，以及相关的MEG解剖数据、3T结构磁共振成像(MRI)数据和空室记录。MEG扫描时长为5-6min，采样率(SR)为2035Hz，抗混叠低通滤波频率为400Hz。

数据处理流程

使用基于MATLAB(The MathWorks Inc. 2022)的Brainstorm软件处理MEG扫描数据，并严格遵循Brainstorm的HCP-MEG教程进行处理。首先利用MEG解剖数据将MEG记录与受试者的结构MRI进行配准。对静息态和空室记录依次施加陷波滤波器(60、120、180、240和300Hz)和高通滤波器(0.3Hz)，分别滤除电源干扰和慢波/直流偏移伪影。随后通过目视检查各受试者记录数据及其通道功率谱密度，剔除不良通道和时间段。利用心电图(ECG)和眼电图(EOG)记录识别心跳和眼动伪影，继而采用信号空间投影(SSPs)方法去除相关伪影。

基于fsLR4k网格定义的8000个源点活动，通过传感器层面记录进行估计。该过程包括：首先采用重叠球模型和垂直于皮质表面的约束偶极子计算头部模型，并根据空室记录估计噪声协方差；随后使用Brainstorm中的动态统计参数映射(dSPM)方法估计源层面活动。最后采用Schaefer-Yeo 7网络100分区图谱对源信号进行分区，其中每个分区的活动是通过计算该分区内所有源活动的主成分来得到的。

利用感兴趣区域(ROIs)的质心坐标定义所有可能区域对之间的欧氏距离(ED)矩阵。区域间距离根据公式δij=EDij/v转换为时间延迟(δ)，其中v为神经信号传导速度(CV)，设定为10m/s。使用Python中的OHBA软件库(OSL)对分区的源定位静息态MEG记录进行源泄漏效应校正。

区域间通信估计

采用两种方法从源层面分区时间序列中推断区域间通信：基于传统FC的方法和事件标记窗口通信(EWC)方法。本研究在FC度量之后添加“-Full”后缀来表示传统的估计技术，使用“-EWC”来表示本研究方法。在估计前，先将源层面分区时间序列划分为10s的时段，剔除坏段，随后对各时段进行独立处理。

“-EWC”方法

本研究开发EWC旨在量化神经元素活动时间序列中离散信号事件或扰动后的功能交互作用。EWC协议包含三个核心要素：(1)事件识别；(2)基于传导延迟的时间排序；(3)事件标记窗口/子样本间的FC估计(图1)。事件识别是指识别活动显著偏离典型行为的特征或实例，具体实现方式是对时段内数据进行z分数标准化：这有助于识别区域活动相对于各自平均行为的离散程度。对于所有区域，将|z|＞3的时间点标记为“超阈值事件”(为简洁起见，下文统称为“事件”)。较高的z分数阈值可筛选出显著强烈的偏差或扰动，本研究认为这些扰动相比于小幅波动更可能对大脑其他区域活动产生可测量的下游效应。有限的信号传导时间还意味着这些下游效应会出现延迟。本研究将区域间通信的估计聚焦于时间序列的特定片段上——称为“通信窗口”，持续时长为1s，源区域从事件时间点开始，目标区域则根据区域间距离换算的时间延迟确定起始时间点(基于区域间距离的时间排序)。一旦确定时间排序的通信窗口，信号事件的下游效应将通过源窗口和目标窗口内活动时间序列的统计关联来量化。EWC具有灵活性，这种关联可使用不同的FC度量方法进行计算，如偏相关(PC)或条件互信息(cMI)。

图1.从区域时间序列估计通信模式的流程示意图。

“-Full”方法

传统方法是在整个时段时间序列上，对所有脑区进行两两之间的FC计算。每个FC估计值均采用α=0.01的阈值进行显著性检验。与EWC方法类似，对所有时段计算两两FC，最终形成贯穿整个扫描时程的N×N×M维矩阵(共M个时段)。

网络模型验证

为了在应用于神经影像数据前评估本方法的优势与局限，本研究在一个具有线性随机模型(LSM)动力学的四节点系统上进行测试。节点i的动态描述为：

其中Cij表示节点i与j之间的连接强度，K为全局耦合参数(K=1)，δij为节点间延迟时间(与欧氏距离成正比)，N(0，1)表示标准正态随机噪声，σ为噪声幅值。其中三个节点以线性链式拓扑结构连接，第四个节点与其他节点隔离，具体的连接矩阵为：C12=C21=C23=C32=1；C13=C31=0；C4i=Ci4=0∀ i∈{1，2，3}。

图2.网络基序上的通信过程。

除了LSM动力学外，泊松过程使两个终端节点以平均频率ν随机产生振幅为0.1的“脉冲”。这导致中心节点的动态是其内部动态与两个终端节点随机输入的混合。虽然泊松尖峰是控制真实信号传递的一种便捷手段，但在高ν值下，由于快速波动可能导致节点平均活动度从“基线”水平虚增。

本研究通过系统改变σ、δ23和ν参数来测试通信协议。对于每组参数值，使用欧拉积分器以2035Hz采样率模拟205s的动力学过程。时段持续时间和通信窗口持续时间分别设为10s和1s。每组模拟进行20个试次。

在所有四个节点的活动时间序列上，分别采用传递熵(TE)和双变量格兰杰因果(GC)的Full方法(TE-Full和GC-Full)，以及条件互信息(cMI)和偏相关(PC)的EWC方法(cMI-EWC和PC-EWC)推断通信。通过判断各方法区分真实连接与虚假连接(相对于真实情况)的能力，即比较1→2和3→2(真实)与4→2(虚假)的估计差异，以及正向与反向连接(如1→2 vs 2→1)的差异，来评估其解析准确通信模式的能力。

计算时间基准测试

为了比较不同协议随网络规模变化的计算时间，对200s时长的时间序列数据分别采用TE-Full、MI-Full、cMI-EWC和PC-EWC推断活动传播。为了确保估计的计算时间的真实性，所用数据为随机抽样的源定位MEG记录。针对网络规模N(共进行10个试次)重复以下步骤：

1、随机选取N个脑区，将其对应区域时间序列截取为200s(采样频率2035Hz)。

2、设置区域间延迟为0.015ms。

3、对所有脑区对采用TE-Full、MI-Full、cMI-EWC和PC-EWC估计活动传播。

4、重新随机选择一个脑区重复试次。

偏相关

偏相关(PC)量化了两个随机变量在排除第三变量影响后的线性关联。在EWC方法中，将第三变量设置为目标区在事件时点前单个窗口时长的历史活动。使用MATLAB的“partialcorr”函数计算PC，显著性水平设为

。偏相关是一个取值范围为-1到1的连续函数，分别表示完全负线性相关至完全正线性相关。

条件互信息

条件互信息(cMI)是一种信息论度量，它量化了在给定外源附加脑区活动(表现为条件变量)的背景下，两个随机变量(代表区域活动)的共享信息量。它是广泛使用的互信息(MI)的一种变体。给定代表两个脑区活动的随机变量X和Y，以及代表第三个脑区(或一组脑区)活动的随机变量Z，cMI通过熵值计算如下：

其中H(X|Z)和H(Y|Z)为变量X和Y的条件熵，H(X,Y|Z)为联合熵。cMI是一个对称测度，取值范围从0(两变量独立)到∞(两变量完全相同)。

格兰杰因果

格兰杰因果(GC)是一种统计预测性度量指标，它量化了一个变量对另一个变量的预测能力。基于向量自回归(VAR)过程定义的双变量GC，通过比较单变量VAR(p)模型与包含另一区域Y历史活动(“因”)的双变量VAR(p)模型对区域X活动的预测误差差异来度量。若引入Y的历史活动能比仅用X自身历史更好地预测X活动(即双变量VAR(p)模型的预测误差低于单变量模型)，则称Y“格兰杰引起”X。为了获得准确有效的GC估计，选择VAR(p)过程的最佳滞后阶数p至关重要。对于实证数据，通过逐步增加滞后阶数拟合模型，估算Akaike信息准则(AIC)等指标，选择对应最小AIC的模型阶数。

对于两个随机且广义平稳的变量X和Y，将 X 的值回归到单变量和双变量的自回归VAR(p)模型上，分别称为约束模型和无约束模型，其回归形式如下：

其中p为模型阶数，X(p)t−1表示X过去值的向量{Xt−1，Xt−2，…，Xt−p}，⨁表示向量拼接，A和A′为VAR模型系数，εt和ε′t为预测误差。X和Y之间的双变量GC定义为：

采用多元格兰杰因果(MVGC)工具箱中高效的状态空间方法估计双变量GC。基于AIC，网络模型的最佳模型阶数选为1，而实证静息态MEG记录的AIC随模型阶数增加单调递减并在50阶左右趋于稳定，故选择更简约的40阶模型(对应约19ms延迟)。

传递熵

传递熵(TE)是衡量两个时间序列信息定向传递的指标。具体而言，对于两个随机变量X和Y，它量化在已知Y自身历史的条件下，X的历史活动能多大程度降低Y(目标)未来的不确定性。从技术上讲，它是cMI的特例，其共享信息量是在相对时移的随机变量间估计的：

其中X(l)t={Xt−l−1，…，Xt−1，Xt}，δ为因果延迟。TE存在多种变体，本研究采用其最简单的形式——双变量TE。多元TE在捕捉真实不对称性方面通常优于其他技术，但需要更大量的数据和计算资源。

如前所述，TE是cMI的特例，因此本研究在事件标记通信窗口内实施的“信号源与目标区的cMI(以目标区历史为条件)”与估计脑区之间的TE密切相关。这体现在二者的公式中：

其中e表示事件时间点，win表示窗口长度。

所有信息论量均使用Java信息动力学工具包(JIDT)的高斯估计器计算。高斯估计器允许基于χ2分布对每个估计值进行显著性检验。除了将TE估计器中的源-目标延迟从默认值1改为基于欧氏距离的区域间延迟外，所有估计器均采用默认参数。所有信息论量(包括EWC)均以纳特(nats)为单位。

结果

本文的分析可分为三个主要部分。首先，在由简单连接模式生成的模拟时间序列中测试了本方法，其中节点活动由线性随机模型(LSM)动力学控制(图2A)。LSM简单的动力学特征使我们能够预先设定节点间的真实信号传递路径。本研究评估了该方法在识别这些通信模式方面的能力，并与先前提出的FC测量方法进行了比较。其次，本研究估算了EWC协议的计算成本，并以常用FC估计方法为基准进行对比。第三，将本方法应用于源定位的MEG记录数据，并且评估了EWC推断的活动传播(FC)与双变量传递熵之间的一致性程度。

网络基序上的非对称信号传递

本研究在具有LSM动力学和泊松尖峰的四节点基序(图2A)中测试EWC协议，以模拟定向信号传递的瞬时事件。EWC采用两种窗内FC测量指标实现——偏相关(PC-EWC)和条件互信息(cMI-EWC)。需要注意的是，传统上这些度量都是对称性指标，无法解析功能交互的方向性。本研究以基于完整时间序列估计的传递熵(TE-Full)为基准，比较这些测量方法在噪声水平递增的模拟中还原真实信号传递模式的能力。为了便于不同测量方法的比较，本研究计算了真实连接(存在：节点1、2与3之间)与虚假连接(不存在：节点4与2之间)以及正向(1→2)与反向(2→1)连接之间FC的标准化对比。标准化对比为零表示估计值之间无显著差异。

本研究考察了方法性能在不同系统噪声水平与信号源区域脉冲振幅比率变化下的表现。信号源(节点1、3和4)与目标节点(节点2)之间的延迟设为15ms，信号源以平均频率为0.2Hz的泊松过程产生脉冲(图2A)。根据实验设计，真实连接存在于节点1与2、3与2之间，而节点4与其他节点没有任何传入或传出的连接。此外，节点1和3的尖峰特征表明它们对节点2存在非对称影响。研究结果发现，通过估算TE-Full，能够很好地反映真实情况，且节点1和3对2的相互作用较为强烈(图2C)。重要的是，相应的对比度分布表明，节点4(孤立)与节点2之间的相互作用可以与真实连接明确区分，证明该方法具有良好的特异性(图2F)。当使用cMI-EWC测量FC时，本研究发现区域4与2之间出现了假阳性相互作用。然而，这种假阳性相互作用明显弱于节点1到2，以及3到2之间的真实连接，尤其是在低噪声水平下(图2D)。有趣的是，PC-EWC的表现优于cMI-EWC。与TE-Full类似，PC-EWC同样能够有效地反映真实连接情况(图2E)。相互作用的对比度显示，即使在高噪声水平下，真实连接仍可与虚假连接区分开来(图2H)，尽管区分度不如TE-Full。值得注意的是，虽然cMI-EWC和PC-EWC都是无定向/对称的FC测量方法，但仍然估算出了不对称的相互作用，这种不对称性源于EWC实现过程中时间顺序的影响。

为进一步评估EWC在解析真实网络通信结构中的性能，本研究采用另一种成熟的定向FC方法，即基于整个时间序列估计的双变量格兰杰因果(GC-Full)。与TE-Full类似，GC-Full也能准确地反映所构建系统的真实相互作用。

此外，本研究还测试了当两个目标区域与共同信号源距离相等时，是否会导致它们之间产生虚假功能连接(闭合三角问题)。为此，本研究将区域2设为脉冲信号源并估计PC-EWC。结果表明，本研究方法正确识别了从节点2到1和3的通信，而反向通信或目标节点间的通信可忽略不计。本研究还测试了EWC在不同传导延迟和信号源平均放电频率下的表现。

综上所述，本研究的对比分析表明，EWC能够在不同噪声水平、延迟和放电频率下准确地捕捉方向性相互作用。需要强调的是，PC和cMI作为无定向/对称的测量方法，仅因EWC协议实现的时序排序才使得这些方法能够捕捉方向性的信号传递。尽管取得了这些积极结果，但基于标准化对比度的观察，本研究发现TE-Full在区分信号传递模式的存在与否方面仍然最为有效。因此，在后续分析中，本研究继续使用TE-Full作为定向FC的代表性基准。

EWC协议的计算可行性

对于任何分析流程(尤其是在探索性研究应用中)，计算可行性都是一个重要考量因素。EWC协议主要包括两个计算步骤，首先是事件识别，然后是通过估计统计依赖性(FC)来推断活动传播。由于一些FC估计方法(特别是涉及非参数信息论函数的方法)本身计算量较大，因此有必要评估新增组件对计算的影响。在这里，本研究比较了EWC协议与TE-Full和MI-Full、以及cMI-EWC和PC-EWC在计算可行性上的差异。使用这些方法估计包含递增节点数量的功能/活动传播网络，节点时间序列采样自源定位的MEG记录。

正如预期所示，所有方法的计算时间均随网络规模增加而增长。有趣的是，尽管EWC增加了计算量(事件识别和针对每个事件的多次FC估计)，但其计算强度仍低于TE-Full。具体而言，在选定参数下，PC-EWC耗时比TE-Full少约75%，cMI-EWC比TE-Full少约45%。虽然方法选择最终需符合具体研究问题，但图3中所示的分析和可视化仍然有助于在静息态分析中根据已建立的方法对EWC协议的计算成本进行基准测试。

图3.网络推断耗时与网络规模的关系。

基于MEG记录推断全脑活动传播模式

在对EWC协议的计算效率和活动传播图推断的准确性进行与传统FC估计方法的对比测试之后，本研究最终测试了从实证神经影像数据中得到的EWC图与TE图之间的一致性。具体而言，本研究使用最小预处理后的30名HCP受试者的静息态源定位MEG数据，计算了PC-EWC和TE-Full。时间顺序基于与区域间欧氏距离成比例的延迟来确定。

本研究计算了每位受试者PC-EWC与TE-Full矩阵之间的相关性，并发现这两种测量方法得出的功能连接估计高度相关(图4B)。受试者间的中位相关系数为r≈0.82(p＜0.0001)(图4C顶部)。值得注意的是，EWC协议的子采样策略显著提高了活动传播推断的速度(图4C底部)。除了在受试者层面上具有高度一致性外，EWC方法还能通过估计不同尺度(如事件、时段等)的非对称程度与估计方差来研究信号传递方向和强度在扫描过程中的动态变化。

图4.EWC在源定位MEG记录中的应用。

总的来说，本研究通过PC-EWC展示了EWC在静息态MEG记录中捕捉活动传播模式的实际应用。研究结果表明，根据EWC协议通过FC推断的活动传播与传统TE方法具有良好一致性，并且计算更高效。这表明EWC是一种可行的技术，特别适用于研究大脑活动或内源性扰动在多个功能性区域间的传递过程。

结论

理解大脑如何传递信息是神经科学领域的一个开放性问题，对于揭示感知与认知的基本机制具有重要意义。本研究聚焦于解决该问题的第一步，即从神经影像数据中可靠地推断通信过程。为此，本研究提出了一种名为“事件标记窗口通信”(EWC)的时间序列分析技术，旨在捕捉由区域间通信驱动的动态统计关系。该方法通过功能连接(FC)度量来捕获离散信号事件或内源性扰动在不同区域间的传播情况。本研究采用不同的FC度量评估了该方法在计算机模拟中的优势与局限，并与现有的方向关系推断方法进行了比较。随后将EWC应用于源定位MEG数据，发现使用EWC获得的受试者水平FC模式与传递熵(TE)、双变量格兰杰因果等方法具有显著的相关性。值得注意的是，EWC具有较高的计算效率(当使用PC时，其每10s时段计算速度约为常规方法的2.5倍)，并且能够捕捉扫描过程中非对称关系的动态变化。综上所述，本研究提出了一种基于活动传播事件的约束FC估计新方法，能够有效推断神经通信模式，该方法为研究大规模神经网络中的信号传递机制开辟了新的途径。

参考文献：Madan Mohan, V., Varley, T.F., Cash, R.F.H., Seguin, C. and Zalesky, A. (2025), Event-Marked Windowed Communication: Inferring Activity Propagation from Neural Time Series. Hum Brain Mapp, 46: e70223. https://doi.org/10.1002/hbm.70223

小伙伴们关注茗创科技，将第一时间收到精彩内容推送哦～