（精读）揭示多元时间序列的高阶组织

原标题《Unveiling the higher-order organization of multivariate time series》

低阶依赖关系： 通过马尔科夫链的模型或分解机FM的模型建模分析；
高阶依赖关系： （跨多个用户项交互的复杂多级级联依赖关系）
（1）基于马尔科夫链方法建模；（2）基于RNN方法；
局限性：模型参数的数量 随阶数呈指数增长，高阶马尔科夫链模型可能涉及的历史状态非常有限，而RNN中采用的过强假设限制了RNN在序列中灵活的高阶应用。
协同波动： 评估矩到矩的协同波动模式（与窗口法相比），量化了各区域在每个时间点上的一致波动幅度。
涌现： 在一定的组织层次上出现的新特性，而涌现特性来自于系统各组成部分的相互作用。

摘要

时间序列分析已经被证明是刻画生物学、神经科学和经济学中的一些现象以及理解它们的一些潜在动力学特征的有力方法。尽管针对多元时间序列的分析已经提出了很多方法，但大多数方法都忽略了非成对相互作用对涌现动力学的影响。在这里，我们提出了一个新的框架来刻画多元时间序列中高阶依赖关系的时间演化。利用网络分析和拓扑学，我们表明，与传统的基于成对统计的工具不同，我们的框架鲁棒地区分了耦合混沌映射的各种时空状态，包括混沌动力学阶段和各种类型的同步。因此，使用模拟动力学过程中的高阶共波动模式作为指导，我们强调并量化了大脑功能活动、金融市场和流行病数据中的高阶模式特征。
总体而言，我们的方法揭示了多元时间序列的高阶组织，从而能够更好地刻画现实世界数据所固有的动态组依赖关系。

正文

来自许多不同复杂系统的丰富的、在时间上可分辨的数据越来越多，这使得详细研究它们的行为和内部机制成为可能。这类系统的例子包括流行病、社会传染、金融、大脑和生物信号。所有这些系统都是由大量的基本单元组成的，它们以异质的方式相互作用，在几乎所有的情况下，都在宏观层面上显示出涌现的特性。由于构成此类系统的相互作用模式至关重要，因此复杂网络成为研究此类系统的结构和动力学的强大框架[1,2]，并有助于表征几种现实世界现象，包括疾病传播[3]，同步[4]，扩散[5]和意见形成[6]。

尽管被广泛认为是许多现实世界复杂系统的参考模型[7-9]，但网络描述的是一个时刻两个单元(或节点)之间的相互作用。然而，这与社会系统[10]、神经科学[11-13]、生态学[14]和生物学[15]中不断增长的群体相互作用的经验证据相冲突。

在上述所有情况下，连接和关系不仅发生在节点对之间，还发生在节点组的集体行动中。通过考虑更精细模型中的高阶(群体)相互作用，如超图和单纯复形[16,17]，最近的几项研究表明，高阶相互作用的存在可以对相互作用系统的动力学产生实质的影响[18,19]，包括改变同步[20]和扩散[21,22]的性质，以及社会[ 23-25 ]和进化过程[ 26 ]中的新集体动力学。

然而，直接测量成对或群体的相互作用来或来约束这些相互作用的这些高阶模型的方法很少。因此，在系统的时空活动模式编码了关于潜在交互的信息的前提下，通常必须依赖于从节点活动的时间序列中提取的间接数据。事实上，在支持多种运动和认知功能的大脑神经元活动中可以观察到这些复杂模式的例子[27,28]，在金融市场中，部分同步模式通常反映金融压力时期[29,30]，也在生物物种的共同进化中[31{33]也有体现。

虽然成对相互作用的推断有很长的历史[ 34 ]，但研究人员最近才迈出重建或过滤高阶相互作用的第一步[ 35-38 ]。特别地，仅依赖成对统计量的方法原则上可能是不充分的，因为重要的信息只能存在于联合概率分布中（考虑两个变量），而不能存在于成对边缘（只考虑一个变量）中，因此无法识别高阶行为[ 39 ]。目前，多变量时间序列中编码的信息在多大程度上源于独立的个体实体，或者更确切地说，源于它们的群体相互作用，仍然是不清楚的。这个问题的一个明显例子是两个大脑区域[ 40、41 ]之间的常规"功能连接"：无论两个区域的活动是成对达到峰值，还是作为一个更大的功能连贯区域的一部分，都可以绘制成对连接。

现有解决这一问题的建议大多局限于描述这种高阶交互的时间性或复杂性，只有少数例外[ 42 ]。例如：
信息论方法通过量化三个或更多相互作用变量组的内在统计协同性和冗余性来表征多元时间序列中的高阶依赖性[ 43-46 ]。
网络科学和随机矩阵理论相结合的一种方法[ 47 ]也被证明适用于揭示相关矩阵的介观组织。
然而，这些方法虽然非常强大，但却难以捕捉到关于系统动态的信息，因为它们需要随着时间的推移进行整合。
例外来自于一个引入了谱分解的工作，该工作将变量组的统计协同和冗余分解为不同的频段，这允许在时间上局部分析高阶依赖性[ 42 ]。
相比之下，来自网络神经科学和信号处理的工具很容易处理多元时间序列的动力学[ 48-55 ]，但只关注成对统计量，忽略了高阶交互作用的影响。因此，仍然缺乏一种原则性的方法来量化节点组的瞬时动态，并可能推断其高阶表示。

在这里，我们提出了一个新的框架来表征信号在所有相互作用阶数(对、三角形等)的瞬时协同波动模式，并研究这种协同波动的全局拓扑结构。我们通过连接时间序列分析、复杂网络理论和拓扑数据分析来实现[ 56 ]。我们首先通过探索时空混沌经典模型所展示的丰富的高维动力学来验证该框架。特别地，我们证明，与传统的基于成对统计的时间序列分析工具[ 57-59 ]不同，高阶测度能够在个体框架水平上揭示不同时空机制的细微差别。然后，我们使用从这些合成模型中获得的见解作为Rosetta stone来解释从时间序列重构的高阶结构，涉及三个不同的真实世界案例研究：静息态大脑活动(通过fMRI数据测量)，股票期权价格和美国各种疾病的流行病学发病率。在所有情况下，我们揭示了额外的丰富的高阶信息，而这些信息在节点和二元网络层面没有被捕获，并突出了不同的拓扑动力学机制，这反过来又会导致系统状态的瞬间分类。

最后，我们展示了推断的动态高阶结构如何提供系统空间配置的瞬时拓扑快照，这些快照可以用作在真实数据集上进一步任务的输入，包括检测大脑区域的局部整合，探索金融危机时期，或从空间传播模式中分类疾病类型。

Results

多元时间序列中高阶结构的拓扑标记

单纯复形非常适合作为描述成对相互作用和高阶相互作用共存的模型框架[ 16 ]。

单纯形（Simplex）：也可称为单形，一个d维单纯形（也表示为d-单形）是由一组(d+1)相互作用的节点形成的。它描述了节点间的多体交互，并允许单形的拓扑和几何解释。例如，节点是0-单纯形，链接是1-单纯形，三角形是2-单纯形，四面体是3-单纯形，以此类推。

单纯复形（Simplicial complex）：又称复形，是一个由单纯形组成的集合，每个单纯形的面都是此集合中的元素，其维度为该集合中最大维度单纯形的维度。

通过单纯复形很容易找出三个元素之间的群相互作用，可以表示为2-单纯形（三个节点），以及节点之间三个成对的相互作用，表示成1-单形（两个节点）。两两相互作用相对于高阶相互作用的相对重要性可以用单纯形的权重来表示，从而形成所谓的加权单纯复形。

我们依靠加权单纯复形表示来描述多个时间序列之间的高阶依赖关系。

算法：
1、对 N 个原始序列进行 z-score；
2、计算所有 k 阶模式的 z-score 的时间序列的元素积，用来表示一个含 k+1 个节点的群体的瞬时协同波动的涨幅；
3、编码 k 阶共波动的时间序列的结果集（即1阶共波动情况下的所谓的边缘时间序列）进一步跨时间进行 z-score，使之在 k 阶上具有可比性；同时选择给产生的权重分配符号，以便在 k 阶序列中区分：（1）完全和谐相互作用（全正或全负）（2）不和谐相互作用（正、负混合）（当 k ≥ 2 时，可能会出现相似的波动，但其机制不同）；
4、对每个时间 t ，将所有瞬时 k 阶共波动浓缩成一个单一的数学对象，即加权单纯复形；
5、对每个时刻 t ，构建一个滤波器 F ( k t ) F(k^{t}) F(kt) ，即按照权重对所有 k 阶波动进行排序，得到一个单纯复形的序列 ∅ = S 0 ⊂ S 1 ⊂ ⋯ ⊂ S l ⊂ ⋯ ⊂ S n ⊂ K t \emptyset=S_0\subset S_1\subset\cdots\subset S_l\subset\cdots\subset S_n\subset K^t ∅=S0​⊂S1​⊂⋯⊂Sl​⊂⋯⊂Sn​⊂Kt 。

Methods

以边为中心：
（1）先对 x i , x j x_i,x_j xi​,xj​ 进行 z − s c o r e z-score z−score;
（2）将两个变量的 z − s c o r e z-score z−score相乘，看共同波动的幅度。
推广到高阶（三角形、四面体 ⋯ \cdots ⋯)
（1）标准化： z i = x i − μ [ x i ] σ [ x i ] z_i=\frac{x_i-\mu[x_i]}{\sigma[x_i]} zi​=σ[xi​]xi​−μ[xi​]​
（2） k + 1 k+1 k+1阶： ξ 0 … k ( t ) = ∏ p = 0 k z p ( t ) − μ [ ∏ p = 0 k z p ] σ [ ∏ p = 0 k z p ] \xi_{0\ldots k}(t)=\frac{\prod_{p=0}^{k}z_p(t)-\mu[\prod_{p=0}^{k}z_p]}{\sigma[\prod_{p=0}^{k}z_p]} ξ0…k​(t)=σ[∏p=0k​zp​]∏p=0k​zp​(t)−μ[∏p=0k​zp​]​
（3）为了区分 k k k阶乘积中的一致相互作用和不一致相互作用，给值赋正负号： s i g n [ ξ 0 … k ( t ) ] = − 1 s g n [ ( k + 1 ) − ∣ ∑ 0 k s g n [ z i ( t ) ∣ ] sign[\xi_{0\ldots k}(t)]=-1^{sgn[(k+1)-\mid\sum_{0}^{k}sgn[z_i(t)\mid]} sign[ξ0…k​(t)]=−1sgn[(k+1)−∣∑0k​sgn[zi​(t)∣]
（4） k k k阶权重 w 0 … k ( t ) = s i g n [ ξ 0 … k ( t ) ] ⋅ ∣ ξ o … k ( t ) ∣ w_{0\ldots k}(t)=sign[\xi_{0\ldots k}(t)]\cdot\mid\xi_{o\ldots k}(t)\mid w0…k​(t)=sign[ξ0…k​(t)]⋅∣ξo…k​(t)∣

为了用本文的方法与基于边缘的方法进行比较，使用边缘时序的均方根代替边缘时序共波动情况， e i j = z i z j e_{ij}=z_iz_j eij​=zi​zj​
R S S ( t ) = ∑ i , j > i e i j ( t ) 2 RSS(t)=\sqrt{\sum_{i,j>i}e_{ij}(t)^2} RSS(t)=i,j>i∑​eij​(t)2 ​

过滤以自上而下的方式进行，从较大的权重到较小的权重—在持续同调[ 61-63 ] —的精神下，这样当k阶单形逐渐包含时，拓扑洞开始出现在 F F F的单形复形中，然后潜在地关闭(即从更相干的模式下降到更不相干的模式)。然而，为了在过滤的每一步保持有效的单纯复形，只能包含满足单纯形闭包条件的k阶单纯形。也就是说，在上一步中子面已经包含在单纯复形中的单形。为了保持这一性质，每当我们添加一个不满足这一要求的单形(例如,一个三角形在它的边之前进入复形)时，我们将其视为简单违背，并将其排除在过滤之外。值得注意的是，这些违反单纯形的结构可以被认为是超相干结构，因为它们的共波动比其子成分(具体见图1d和方法)的共波动要强。更进一步地，在本文中我们给出了当 k = 2 k = 2 k=2 时的结果，因此我们只考虑到三角形以内的简单情况。尽管如此，我们的框架自然地推广到更高阶的(即k≥3 )。

Methods of hyper coherence and hyper complexity

思想： 构造一系列的单纯复形，以越来越高的精度逼近原始的加权单纯复形， ∅ = S 0 ⊂ S 1 ⊂ ⋯ ⊂ S l ⊂ ⋯ ⊂ S n ⊂ K t \emptyset=S_0\subset S_1\subset\cdots\subset S_l\subset\cdots\subset S_n\subset K^t ∅=S0​⊂S1​⊂⋯⊂Sl​⊂⋯⊂Sn​⊂Kt
（1）以递减的顺序对链式和三角形权重进行排序；
（2）在每一步 l l l中去掉所有不满足单纯闭包条件的三角形，即 w i j < w i j k w_{ij}<w_{ijk} wij​<wijk​,也就是说，该三角形里有比整个三角形权重小的边，那个这个三角形就是违规三角形，去掉该三角形，并在违规行为列表中插入相应的权重 Δ v = { ( i , j , k ) ， w i j k } \Delta_v=\lbrace(i,j,k)，w_{ijk}\rbrace Δv​={(i,j,k)，wijk​}，其余权重比该权重大的链、三角形都属于单纯复形。

过滤与持久同调性
1D-循环是没有任何分支的闭合环路
1D-循环的持久发生器：第一个同源群 H 1 H_1 H1​的提供了关于何时何地出现更高同步区域的见解。
（a）二维简单过滤；
（b）长条描述了不同拓扑特征在多个尺度上的生命周期，每个长条对应于一个特定的拓扑特征，可以根据时间戳识别，“时间戳”记录了该特征在过滤过程中的开始与结束。绿色： H 0 H_0 H0​各种连接组件的持久性，他们逐渐合并，直到只有一个存活；蓝色： H 1 H_1 H1​描述一维循环的生命周期；
（c）持久性图和条形图等效，例如，若只关注 H 1 H_1 H1​，则每个一维循环的二维图中由坐标 ( w b , w d ) (w_b,w_d) (wb​,wd​)进一步区分，一维循环的性质取决于相应的时间戳，即一维循环是否在 w ∗ = 0 w^*=0 w∗=0之前还是之后创建与关闭，它反映了 F D 、 C T 、 F C FD、CT、FC FD、CT、FC结构。（完全相干FC： 创建和关闭一维循环边和三角形的 w > w ∗ w>w^* w>w∗；完全退相干FD： 创建和关闭一维循环边和三角形的 w < w ∗ w<w^* w<w∗；相干迁跃CT： 创建一维循环的 w > w ∗ w>w^* w>w∗，关闭一维循环的 w < w ∗ w<w^* w<w∗。）
（d）超复杂性指标定义为 H 1 H_1 H1​的持久图与空持久图的 W a s s e r e s t e i n Wasserestein Wasserestein距离，可以根据一维循环性质分解成三部分，并绘制成三角形图。

超相干： 违反相干三角形 ÷ 相干三角形
平均违反指标：违反三角形上平均的边总数（在进入三角形前有几条边在单纯形中）
持久同调： 持久性 π g = d g − b g \pi_g=d_g-b_g πg​=dg​−bg​ 死亡的时间-出生的时间
H 1 H_1 H1​持久同源群描述的是一维循环在过滤过程中存活的时间，作为最底层空间“复杂性”代理， H k H_k Hk​同调生成的持久性之和等于拓扑空间到平凡空间的距离（不包含 k + 1 k+1 k+1维空间）。
超复杂性： H 1 H_1 H1​的持久图——空持久图
同源支架（Homological scaffold）： 更好的描述持久性图中存在的拓扑特征。
该对象是一个由生成元 g i g_i gi​所对应的所有循环路径组成的加权网络，其中生成元 g i g_i gi​的权重由它们的持久性 π g i π_{ gi} πgi​决定。换句话说，如果一条边e属于多个一维圈 g 0 ， g 1 ， … ， g s g_0，g_1，\ldots，g_s g0​，g1​，…，gs​，则它的权重 w ˉ e π \bar{w}^ π _e wˉeπ​ 定义为生成元的持续度之和，即： w ˉ e π = ∑ g i ∣ e ∈ g i π g i \bar{w}^ π _e=\sum_{g_i\mid e\in g_i}{\pi_{gi}} wˉeπ​=gi​∣e∈gi​∑​πgi​同调支架提供的信息允许我们破译不同的链接在系统的同调性质方面所起的作用。一个链路 e e e的总持续性 w ˉ e π \bar{w}^ π _e wˉeπ​ 意味着该链路在相干和非相干协涨落空间中充当局部强桥[ 11 ]。
低阶投影（lower-order projections）： 最后，为了分析边/节点级别的违规列表 Δ v \Delta_v Δv​所提供的信息，我们采用向下投影的方法。也就是说，对于每条边 ( i , j ) ( i , j) (i,j)，我们分配一个权重 w i j w_{ij} wij​，等于由该边定义的三角形的权重和的平均，即具有权重 w i j ⋅ w_{ij \cdot} wij⋅​的形式 ( i , j , ⋅ ) ( i , j , \cdot) (i,j,⋅)的三角形，其平均值是在由该边定义的三角形的数量 n i j ⋅ n_{ij\cdot} nij⋅​上计算的。类似地，我们定义节点 i i i 的节点强度 w i w_i wi​为与节点 i i i 相连的三角形的权重之和的平均值。在同调支架的情况下，由于它是一个加权网络，节点 i i i的节点强度 w i ˉ \bar{w_i} wi​ˉ​用经典的[ 8、116 ]方法定义为与节点 i i i 相连的边的权重之和。

总之，对于每个时刻 t t t ，该框架会产生两种结果：
1、有简单的闭合条件诱导出一个违背三角形列表，即 Δ v = { ( i , j , k ) , w i j k } \Delta_v=\{ (i,j,k),w_{ijk}\} Δv​={(i,j,k),wijk​}，这些是2—单纯形(三角形)，它们的权重至少比其对应的1 —单纯形(边)的权重大。从直觉上看，这些三角形反映的是高阶状态，不能仅仅通过成对的协同波动来捕捉。然后，我们定义了超相干指标，即违反相干三角形(即权重大于零的违规三角形)在所有可能的相干三角形(即权重大于零的三角形)中所占的比例 。
2、单纯形过滤 F F F是嵌入的单纯形复合体序列- -按照相干模式排序- -以空复合体为起点，以整个单纯形复合体(见图1d的右面板)为终点。然后，我们通过计算 F F F的持久同调来刻画某些拓扑特征(连通分支, 1 -维圈, 3D -洞等。) [ 63、64 ]的持久性质。在这里，我们重点研究过滤中的1D循环，即第一同调群H1的持续生成元，它提供了关于何时何地出现更高的同步区域的见解。持久同调的经典输出是一个条形码(或者等价地,一个持久性图)，它是一个压缩的总结，描述了1D循环在 F F F (见图S1)上存活的时间。我们以此为对象，将超复杂度指标定义为H1的持久图与空持久图之间的Wasserstein距离[ 65 ]，对应于一个H1同调(详见Methods)的平凡空间。通过这种方式，我们得到了相干和非相干共波动的拓扑复杂性的度量。

全局和局部拓扑标记对不同的动力学机制进行分类

为了深入了解我们的拓扑指标的性能，我们在这里证明了超相干性和超复杂性容易区分由时空混沌的典型模型产生的不同动力学区域。作为一个案例研究，我们考虑扩散耦合映象格子( CMLs ) [ 67 ]，它是定义在离散时间和空间上的具有连续状态变量的高维动力系统。CMLs被广泛用于多个不同领域的复杂时空动态建模，包括生物学[68]和金融[69,70]。特别地，我们考虑一个具有 N N N 个格点的环形格点，并假设每个格点 i i i 的状态 x i x_i xi​ 的系统的动力学演化是两个不同的竞争动力学的结果：内部混沌动力学和第一近邻格点之间的外部扩散耦合动力学。它们的动力学可以表示为 x i ( t + 1 ) = ( 1 − ε ) f [ x i ( t ) ] + ( ε 2 ) f [ x i − 1 ( t ) + f x i + 1 ( t ) ] x_i(t+1)=(1-\varepsilon)f[x_i(t)]+(\frac {\varepsilon} {2})f[x_{i-1}(t)+fx_{i+1}(t)] xi​(t+1)=(1−ε)f[xi​(t)]+(2ε​)f[xi−1​(t)+fxi+1​(t)] 其中 i = 1 , 2 , … , N i=1,2,\ldots,N i=1,2,…,N , ε ∈ [ 0 , 1 ] \varepsilon \in [0,1] ε∈[0,1] 是耦合参数（或耦合强度）， f ( x ) f ( x ) f(x)一般为混沌映射。在所有的模拟中，我们都考虑了Logistic映射，即 f ( x ) = 1 − 1.75 x 2 f ( x ) = 1 - 1.75 x^2 f(x)=1−1.75x2。

现在已经很好地确定了[ 66,71 ]，通过改变混沌映射的耦合强度 ε ε ε 的值，CMLs表现出各种各样的时空模式，包括不同程度的同步和动力学阶段，如

充分发展的湍流( FDT ）：一个具有非相干时空混沌和高维吸引子的阶段)，
模式选择( PS ）：一个有利于随机选择周期吸引子的混沌抑制阶段，反映了准周期行为，
不同形式的时空间歇( STI）：在FDT和PS之间插入的具有低维吸引子的混沌伪相位；
缺陷布朗运动( BMWD）：系统中存在缺陷并于布朗运动类似混沌波动的阶段，
缺陷湍流( DT）：一个产生许多缺陷并湍流地碰撞在一起的阶段)[ 72 ]。

值得注意的是，这个非常丰富的相图的起源来自于由每个单一状态的混沌动力学引起的局部趋向不均匀性和由扩散动力学引起的系统在空间均匀化的全局趋向之间的相互作用[ 72 ]。

figure2：全局和局部高阶指标区分了耦合混沌映射的动力学状态。(a)我们报告了一个 N = 119 N = 119 N=119个节点， T = 1200 T = 1200 T=1200的多元时间序列的超相干性指标的时间演化，该指标是由固定时间长度 L = 240 L = 240 L=240的5个不同的CMLs机制串联得到的。即，从有序到无序， ε = 0.12 ε = 0.12 ε=0.12时的模式选择( PS )， ε = 0.08 ε = 0.08 ε=0.08时的带缺陷的布朗运动( BMWD )， ε = 0.3 ε = 0.3 ε=0.3时的时空间歇性II ( STI )， ε = 0.068 ε = 0.068 ε=0.068时的缺陷湍流( DT )和 ε = 0.05 ε = 0.05 ε=0.05时的充分发展湍流( FDT ) [ 66 ]。(b)值得注意的是，当将违反三角形 Δ v Δ_v Δv​的列表投影为加权图(见向下投影定义的方法)时， 边权分布 P ( w i j ) P ( w_{ij} ) P(wij​)反映了从有序到无序的不同动力学机制的性质和"排序"。（c)我们绘制了超复杂度指标的时间演化图和,(d)由H1的同源生成器构建的同源支架的权重 P ( w ˉ i j ) P (\bar{w}_{ij} ) P(wˉij​)的分布[ 11 ] 。作为比较，我们还在面板 ( a , c ) (a,c) (a,c)中报告了在没有任何约束的多变量时间序列(灰色曲线)重新洗牌时获得的零模型的相同指标。关于更保守的零模型中高阶指标的行为，见SI第S4节。阴影区域和误差条代表100个独立实现中的标准偏差。

在图2中，我们总结了我们的高阶方法应用于这些N = 119节点和T = 1200的合成多元序列时的结果，这些序列是通过连接固定时间长度L = 240的cml的五个不同动态阶段获得的。即，从有序到无序，其中 1 0 5 10^5 105个时间点的瞬态已被去除。在SI图S2中报告了这种多变量时间序列的样本，而我们在SI Section S2.2和SI图S3中研究了z-score的影响。值得注意的是，图2(a)中报道的全局超相干指示器清楚地区分了cml的不同动态相位，同时也保留了有序和无序状态之间的排序。更准确地说，它赋予完全和部分同步状态高值，而相反，混沌或湍流状态表现出较低的超相干值。虽然该指标仅提供全局信息，但可以通过将违规三角形列表 ∆ v ∆_v ∆v​的大小投影为加权图来获得精细信息(见定义向下投影的方法)。同样在这种情况下，事实上，边权分布 P ( w i j ) 反 P(w_{ij})反 P(wij​)反映了不同动力状态的性质和秩(图2b)。周期序列，例如PS，可以转换为峰值良好的分布，类似于泊松分布。而当无序度进入多元时间序列的伪相时，边缘权值分布逐渐改变形状，FDT混沌序列的极限情况收敛于不均匀分布。

在研究超复杂性的时间演化时，也可以得出类似的结论(图2c)。然而，我们发现了一些显著的差异。对于复杂性，最低值分配给周期性模式(例如PS)，因为这些制度需要少量的信息来描述。相反，FDT等混沌状态显示出最高的超复杂性值。虽然这个高阶指标也能够区分 C M L s CMLs CMLs的不同动力学机制，但人们可能会认为，超相干性和超复杂性指标提供了由强负相关(即斯皮尔曼秩相关 ρ ≈ − 0.95 ρ≈−0.95 ρ≈−0.95)所表明的等效信息。我们将在下一节中说明，对于实际的多变量时间序列，这通常是不正确的。

最后，在图2d中，我们报告了持久性同源支架的边权分布，这些图是由H1的持久性同源生成器构建的(详见Methods和参考文献[11])。这些分布量化了共波动中边的拓扑重要性，即它们所属的同源发生器的持久性。值得注意的是，当我们从周期到混沌多变量时间序列时，这些分布也改变了它们的整体形状，在这种情况下，保持有序和无序之间的秩。

z-scores和符号分配对CMLs时序的影响

为了区分k阶乘积中的完全协和群作用(所有正向或负向的波动)和不协和群作用(正负波动的混合)，我们给完全协和群作用分配正符号，给不协和群作用分配负符号。然而，这种分配符号的方式意味着基线分数可能会稍微改变，因此，三角形的一致性可能比边更难实现。

边权分布—三角形权重分布—无符号分配的三角形权重分布

我们分析了边和三角形分布中 z − s c o r e z-score z−score和符号分配对耦合映射晶格动态状态的影响。特别是，三角形的 z − s c o r e z-score z−score和符号分配相对于边的 z − s c o r e z-score z−score和符号分配产生负偏移。这种现象对超相干指示器没有很大影响，因为它被定义为违反相干三角形(即违反权值大于零的三角形)的比例超过所有可能的相干三角形(即权值大于零的三角形)。相比之下，三角形权重分布在边缘上的负偏移量可能会对超复杂性指标产生影响。这是因为该指标是通过考虑边和三角形权重分布的集合来构建的。然而，对超复杂性的影响是贡献从完全相干到其他两种贡献的转变，而不是整个简单过滤的重组。从这个意义上说，这种效应是有偏差的。

figure3。**高阶方法在区分CML状态方面表现较好。**对于五种动力机制，我们报告了(a)超相干，(b) R S S RSS RSS[54]， © S − i n f o r m a t i o n S-information S−information[43,45]和(d) P e r s o n Person Person相关分布的小提琴图。值得注意的是，只有高阶方法能够区分五种动力状态(即 I C C > 0.9 ICC > 0.9 ICC>0.9)。还要注意，静态方法只能在对每个块的位置有先验知识的情况下使用。

这些结果定性地证实了用我们的方法提取的全局和局部拓扑信息都能很好地区分不同的动态状态。我们用类内相关系数(intraclass correlation coefficient, ICC)[73,74]定量评估了高阶指标区分动态机制的能力，类内相关系数是一种常用的统计度量，用于确定不同群体的单元(或评级/分数)之间的一致性。换句话说，ICC描述了同一组中的单元彼此相似的程度，因此协议越强，其ICC值就越高。在图3中，我们报告了几种方法在试图区分CML的五种动态机制时的比较。从理论的角度来看，有趣的是，从图3a中可以看出，五种 C M L s CMLs CMLs体系映射成不同的超相干分布，反映了每个时空体系固有的一些特征。例如，除了定义良好的体外，STI区还呈现出一条长尾，指向较低的超相干值，这捕获了由于动力相位不匹配而引起的罕见的短混沌爆发[66,75]。

我们发现我们的两个高阶测量(在图3中只显示了超相干性)都有很高的 I C C ICC ICC值(即，对于超相干性和超复杂性，分别约为0.95和0.96)。对于一种常用的时间低阶测量， R S S RSS RSS统计量[54]，它解释了所有一阶共波动的峰值幅度，即边缘时间序列(见正式定义的方法)，我们发现ICC值相当小( I C C ≈ 0.57 ICC≈0.57 ICC≈0.57)，如图3b所示，这意味着高阶效应占主导地位。事实上，我们的动态方法与“静态”高阶信息理论方法相当[43,45]，即对于s-information(如图3c所示)，我们发现 I C C ≈ 0.98 ICC≈0.98 ICC≈0.98。然而，这些数量通常是在时间窗口上计算的，而在这里提出的拓扑框架中，有可能具有瞬时信息。最后，正如预期的那样，基于Pearson的相关性(即 I C C ≈ 0 ICC≈0 ICC≈0)，高阶度量(静态和动态)始终优于静态低阶方法，如图3d所示。事实上，显然只有高阶方法才能有效地区分各种时空状态，而低阶统计无法捕捉动态状态之间的细微差异。与其他“静态”高阶和两两方法的详细比较[41,47,54]，参见SI章节S3和SI图S6-S7。

不同零模型的影响

（1）块零模型： C M L s CMLs CMLs每个动态机制对每个区域（240个时间点组成的块）进行洗牌，重组的只是每个机制内部的数据；
（2）相位随机化零模型：利用傅里叶变换将 C M L s CMLs CMLs时间序列变化到频域上，再对相位系数洗牌，再变到时域上。

（1）块零模型中，超相干块变得彼此难以区分，同时保持了共同波动的某些特定状态特征。没有保持有序状态到无序状态。超复杂性指标几乎相同，不能在块零模型中分离；（2）相位随机化零模型时重复相同的分析，该模型保留了CML时间序列的功率谱。

耦合映射格中空模型的比较。(a-d)我们报告了经典CML时间序列的超相干性分布，以及本工作中考虑的三种零模型，即正文中考虑的简单零模型、块零模型和相位随机化零模型。(e-h)对ICC度量的全局和两两分析都表明，只有块零模型能够部分分离CML区域，从而保持了共同波动的某些特定区域特征，同时打破了真实超相干分布中存在的有序-无序梯度。(i-l)在重复分析超复杂性分布的面板(a-d)时，我们发现，对于五个区块，零模型中超复杂性的分布大多是相同的，因此很明显，存在一些仅存在于CML多元时间序列中的固有拓扑属性。(m-p) ICC度量的全局和两两分析提供了一种定量的方法来评估不同零模型中每个制度之间的相似性。

现实世界的复杂系统表现出非平凡的超相干性结构。

作为对真实世界多元时间序列分析的应用实例，我们报告了高阶框架在人类连接组计划( HCP ) [ 76 ]的fMRI信号、纽约证券交易所( NYSE )的金融资产价格以及美国[ 77、78 ]的几种传染病历史数据上的结果。

对于人脑数据，我们考虑HCP 100名无关被试的静息态fMRI信号，采用HCP发布提供的100个脑区[ 79 ]和19个皮层下脑区[ 80 ]进行皮层分割，共计N = 119个感兴趣区域( Regions of Interest，ROI )。对于金融时间序列，我们分析了2000 - 2021年纽约证券交易所部分美国公司N = 119个股票价格的逐日时间演变。最后，对于流行病数据集，我们调查了衣原体、淋病、流感、麻疹、腮腺炎、脊髓灰质炎和百日咳(详见数据集上的方法)在美国州级( N = 50)的每周病例数。

figure4.面向真实世界多元时间序列的高阶指标。( a )小提琴图显示了三个真实世界数据集的超相干性分布，即静息态fMRI数据( N = 119个脑区)，纽约证券交易所119种资产的金融价格，以及美国国家层面的几种传染病的历史数据( N = 50 )。将真实分布与5种CML动力学机制，以及在没有任何约束的独立重排时获得的相应零模型进行比较，合成和真实的多变量时间序列(见SI节S2.3 ,图S5 ; SI部分S4 ,图S8 - S10为更多细节)。注意到所使用的三个真实世界数据集的分布在其轮廓中显示出明显的差异，但在统计上总是与相应的零模型不同。见SI剖面图S5和SI图。( b )相干和非相干共波动模式中与一维循环相关的不同贡献的二维直方图。在这里，三角形中每个点的位置由与一维循环相关的三个不同贡献决定。例如，如果将超复杂度指标拆分为完全相干( FC )、相干转移( CT )和完全退相干( FD )三个相等的贡献，则一个点将位于三角形的中心，而角点的位置保留给那些主要贡献来自FC、CT或FD的点。

在图4a中，我们报告了五个CML动力学区域和三个数据集的超相干分布。为了比较，我们还绘制了独立重排合成和真实世界多元时间序列(更保守的零模型中高阶指标的行为见SI Sections S4和SI Figures S8 - S10)得到的零模型。考察现实系统的超相干性分布时，可以观察到几件事情。首先，这些分布总是与相应的零模型(所有 p < 1 0 − 10 p< 10^{-10} p<10−10 ,采用科尔莫戈罗夫- Smirnov检验)在统计上截然不同，但它们也表现出彼此截然不同的特征。如果我们关注流行病数据，例如，已经可以通过粗略地比较相应的超相干分布来区分疾病。这些分布实际上反映了疾病演化所固有的独特的高阶时空模式。对于金融系统，通过对比，我们得到了双峰分布，反映了危机和稳定的金融时期之间的二分法。也就是说，经济危机的典型特征是增加(超)同步性，而金融稳定时期似乎以更混乱的方式展开。此外，结合CMLs的解释基准，我们发现，在休息期间，人类大脑主要与混沌状态相关，少数与部分同步状态相关，这与静息态大脑动力学的研究一致[ 82-85 ]。

超复杂度分解提供了关于动态机制的详细信息。

为了更好地刻画一维同调生成元在相干和非相干共波动模式的演化，我们将超复杂度指标分解为三种不同的贡献。也就是说，当我们追踪一维周期沿过滤的演化时，我们关注的是仅由完全相干结构创建和闭合的一维周期，即权重大于零的边和三角形，我们称之为全相干( Full Coherence，FC )贡献；由相干结构形成的一维循环被非相干结构(即权重小于零的边和三角形)闭合，我们将其定义为相干跃迁( Coherence Transition，CT )；最后，仅由完全退相干结构产生的1D环，我们称为完全退相干( FD )贡献。显然，通过构造，这三个贡献的总和可以归结为总的超复杂性。我们在SI图S1中展示了一个示例。

在图4b中，我们将三个对超复杂性指标的贡献绘制在一个三角形表示中。在这个空间中，如果所有的一维循环都是由完全相干结构形成和闭合的，那么一点被放置在左下角。同理，右下角对应完全非相干结构的排他性贡献，上角对应由相干跃迁唯一决定的贡献。每当超相干指标分裂为相似的FC和FD贡献时，该点被放置在相应的角点之间，因此其位置反映了贡献的相对重要性。例如，如果将超复杂度指标拆分成三个相等的FC - CT - FD贡献，则一个点将位于三角形的中心。注意到这种分解携带了关于超相干性指标的完全不同的信息，但我们得出了与刚刚提出的结果类似的类比。

事实上，当考察超复杂度指标在合成信号中的不同贡献时，我们发现五个 C M L s CMLs CMLs机制在不同的簇中似乎是分开的。部分同步信号由全相干性和相干跃迁贡献的混合来表征，而混沌信号主要由全相干性和非相干性(见图4b左图)决定。相比较而言，对于静息状态下的人脑，我们发现大部分的状态被定位在混沌和部分同步状态之间。这与考虑超相干性指标时得到的结果一致，该指标提供了不同性质的信息，即它仅基于单纯的违规次数。

高阶拓扑标记的实际应用

到目前为止，我们主要关注了我们的全局高阶指标在合成和真实世界的多元时间序列中的时间演化。下面我们报告在更局部的层面上考虑高阶测度时的一些代表性应用。我们的目标是刻画在静息态脑数据和金融系统中具有最大同步水平的高阶状态。为此，在人脑的背景下，我们分离出了前15 %的相干帧，这些相干帧与更同步的动力学相位有关。在图5a中，我们通过将违反三角形 Δ v \Delta_v Δv​的大小投影到一个节点级别(详见方法和SI图S15 ,用于其他峰百分比的比较)上，报告了一个最具判别性的节点的脑图。这相当于考虑了从违反三角形 Δ v Δ _v Δv​列表中提取的节点强度。换句话说，绝对值最大的区域是属于最连贯的高阶结构的区域。特别地，我们发现具有强调同步协同波动的活动模式主要反映感觉运动区域，它们属于存在于静息态网络中的著名底物之一[ 86 ]。当考虑直方图报告七个典型功能网络的平均相干性时，这一点得到了证实[ 81 ] (也见SI图S16为高阶指标对功能网络之间的影响)。

figure5。高阶测度的投影提供了局部时空信息。从从违规三角形 Δ v \Delta_v Δv​中提取的节点强度可以用来跟踪高阶结构在时间上的重要性。( a )当分离前15 %的相干帧时，得到了参与高阶协同波动的节点的脑图，这些节点是与休息期间更同步的动态相关的节点。我们还展示了报告Yeo七个典型功能网络[ 81 ]中平均相干性的箱线图，即视觉( VIS )，躯体运动( SM )，背侧-注意( DA )，腹侧-注意( VA )，边缘( L )，额顶( FP )和默认模式网络( DMN )。( b )产业部门层面违反三角关系的节点强度 Δ v \Delta_v Δv​的时间演变将危机与金融稳定时期区分开来。( c )当选择前15 %低的超复杂特性，获得的脑图主要包括默认模式网络，这被其高平均节点得分所证实。( d )同调框架节点强度的时间演化提供了某些经济部门衰退的更精细的细节。值得注意的是，从拓扑角度来看，从同调支架中提取的节点强度提供了关于共波动空间中一维循环的信息。( e )最后，通过 t − S N E t - SNE t−SNE非线性降维，将超相干性、超复杂度的三个贡献和平均边违反度(见定义法)作为特征，可以得到流行病爆发历史数据的平面嵌入。作为内插图，我们报告了在选择15 %高相干框架时，美国国家层面违反三角形 Δ v \Delta_v Δv​的节点强度。值得注意的是，疫情的时空演化在不同状态和不同疾病之间存在差异。

在图5b中，我们报告了从违反三角形 Δ v \Delta_v Δv​列表中提取但在工业部门层面聚集的节点强度对金融时间序列的时间演化。最高值捕捉到了金融不稳定的主要时期( 2002年对应市场低迷, 2007 ~ 2008年对应次贷危机导致的大衰退)的开始，其特征是股票价格的同步性增强，这与2002 - 2007年和2013 - 2018年的不同步区间明显区分开来，而这反过来又对应着经济更稳定的时期。

类似的分析可以通过关注H1(具体方法参见)的持续同源生成器构建的同调支架的超复杂度指标和节点强度来产生。特别地，图5c描述了当分离15 %的低-超复杂度帧时得到的大脑图，如前所述，这些帧是与更同步的动态相位相关的帧。在这里，绝对值最高的是与默认模式网络( DMN )相关的网络，DMN被认为是唤醒休息[ 87、88 ]期间最活跃的网络。

相比之下，对于图5d中的金融时间序列，同调框架节点强度的时间演化提供了某些经济部门衰退的精细细节。例如，消费品、基础材料以及石油和天然气是受2007年大衰退影响的主要部门。

最后，通过分析美国流行病爆发的历史数据，表明高阶测度(即超相干性、超复杂度的3个贡献和平均边缘侵犯;定义见方法)的时间演化可以用于不同传染病的分类。特别地，使用重复50次(分类器之间的比较见SI表S2)的10倍交叉验证设置，支持向量机( SVM )分类器报告了较高的准确率水平，约为85 %。为了更直观地表示这一结果，我们在图5e中报告了使用t分布随机邻域嵌入( t-SNE )非线性降维方法获得的流行病爆发历史数据的平面嵌入。值得注意的是，非线性方法，如t - SNE，允许在投影到低维空间后保留原始高维空间中的"局部"结构，这通常是主成分分析( Principal Component Analysis，PCA )或多维尺度分析( Multidimensional Scaling，MDS )等线性方法无法做到的[ 89 ]。在这个空间中，我们观察到不同种类的疾病在很大程度上聚集在一起，这在一定程度上反映了疾病爆发的独特时空演化特征，这些特征被SVM分类器捕获。同时，可以观察到疾病之间的相似性。对于性传播疾病，如淋病和衣原体，这种情况在平面嵌入中大多是重叠的。作为插图，我们还报告了选取15 %高相干框架并考虑违反三角形 Δ v \Delta_v Δv​的节点强度时得到的美国国家层面的地图。我们发现，疾病暴发的时空演化在不同状态和疾病之间存在差异，一定程度上反映了疾病独特的"高阶"特征。

讨论

推断多元时间序列中高阶结构的动力学性质在许多复杂系统中至关重要，从流行病学，到金融，再到生物系统。然而，直接的高阶网络测量往往是无法访问的[ 18 ]。事实上，在许多生物、社会和金融系统中普遍存在的绝大多数复杂的时空活动模式通常是在节点层面上记录的，而不是在边或组的层面上直接测量的。本文介绍的高阶方法提供了第一种从多元时间序列动态重建高阶相互作用的强大和替代方法。

作为一个初始基准，我们首先验证了我们的方法针对那些基本动力学已知的信号。特别地，与各种低阶统计量[ 47、54 ]不同，本工作提出的全局高阶指标能够鲁棒地分离高维耦合混沌映射中的若干动力学相位，这些动力学相位似乎只能通过基于高阶统计量的方法来区分[ 43 ] (见图3和SI图S6、S7)。这为识别高阶行为的高阶方法的必要性提供了进一步的经验证据[ 39 ]。在更局部的层面上，即当把超相干三角形列表投影为加权图时，图的权重分布反映了多元时间序列的全局动态：同步的周期序列转化为尖峰分布，而混沌序列转化为胖尾分布。有了这些理论基础，我们将我们的框架应用于现实世界的多元时间序列，特别是来自100个不相关的人类受试者的静息态fMRI信号，纽约证券交易所的金融资产价格，以及来自7个不同疫情爆发的历史数据。

我们发现，在休息期间，人脑的高阶动力学主要在充分发展的湍流和部分同步之间振荡。这与最近的研究一致，支持人类大脑工作在一个动荡的状态[ 82、90 ]，处于临界性的边缘[ 91 ]，这似乎赋予了显著的信息处理优势[ 92 ]。此外，当在更精细的尺度上分析大脑状态时，我们发现了两个没有表的方面。一方面，最大相干高阶结构反映了感觉运动区域，属于静息态脑网络结构中存在的已知底物之一[ 86 ]。众所周知，这些区域在破译- -在非常快的时间尺度上[ 55 ] - -由于外部环境而不断变化的输入方面发挥了主要作用[ 83 ]。另一方面，在考察其最低点的超复杂度标记时，我们发现节点投影粗略地捕捉到了默认模式网络( Default-Mode Network，DMN )，该网络已知整合了人脑网络[ 87、88 ]中的高、低阶信息。因此，我们提出的两个标记物提供了互补的见解，这些见解不是简单地从摆振倾角方法(见SI断面S5和SI图S12、S13)推导出来的，关于大脑网络是如何随时间分离和整合高阶信息的：超同步的，较少整合的，这是由感觉运动区域的简单侵犯的数量来衡量的；以"恰当"的高阶简化来衡量，在系统内部整合度更高的DMN，其脑动力学主要作用于内向，通过整合低阶和高阶动力学[ 55、93 ]，处于"不稳定的边缘"，处于内部探索的恒定状态[ 84 ]。

相反，在金融时间序列的背景下，我们提供了证据，高阶结构的幅度有效地区分了危机和金融稳定时期，这不能从不同的零模型(见SI节S4.2)中获得。特别地，最大一致高阶结构在重大金融危机对应中强烈涌现，反映了同步共激活模式(即股票价格倾向于向同一方向移动,因此提高了它们的同步性水平)的增加。虽然这在文献[ 30、72、94]中并不是新的，但我们强调，与我们的方法不同，大多数这些方法依赖于在滑动时间窗口上估计的相关矩阵[ 29 ]，因此忽略了在个别帧(例如在高频交易中)水平上可能要捕获的信息。在考察处于最低点的超复杂性指标时，拓扑标记捕捉了关于工业部门在危机期间的不同作用的精细信息，揭示了不同行业在时间上的强烈差异和高度异质性[ 96 ]，表明它们具有识别系统性风险[ 97、98 ]累积的潜力。

最后，在分析美国流行病爆发的历史数据时，我们已经表明，我们的高阶测度的时间演变可以在很大程度上对不同种类的疾病进行分类。特别地，疾病的平面嵌入揭示了基于其独特的时空模式的有趣集群的存在。虽然这个结果本身在不同的学科中是有趣的[ 99-102 ]，我们的高阶标记可能为探索流行病爆发的可预测性[ 77、103、104]提供新的工具，尽管当增加预测长度[ 77、105 ]时，它的局限性仍然存在。

综上所述，这里我们开发了一个新的灵活的工具来提供多元时间序列中高阶结构的逐帧估计(因此规避了滑动窗口方法的局限性[ 106 , 107 ])。我们相信，我们的框架可以有效地用于所有对信号的动力学知之甚少或未知的情况，为进一步在生物学、流体动力学、社会科学或临床神经科学等领域的应用铺平道路。特别地，高阶指标和相应的低阶投影还可以提供拓扑Polaroids，即所研究系统空间构型的瞬时拓扑快照。总的来说，我们的方法表明，与标准方法相比，研究多元时间序列的高阶结构可能提供新的见解，从而更好地刻画现实世界数据固有的组依赖关系。

限制

我们的方法的主要限制之一是时间复杂性。事实上，如果我们考虑到k阶的共波动模式，得到的时间复杂度为 O ( N k ) O(N^k) O(Nk)。此外，在当前阶段，我们的框架不允许调查两个后续时间框架之间的因果关系(即，就所提出的拓扑标记而言，之前的时间点对下一个时间点的影响程度)。还要注意，我们的动态高阶方法，和许多其他现有的两两动态方法一样[48,49,54]，会受到时间序列中有噪声波动的严重影响。然而，这个问题可以通过分析多个时间框架的平均统计数据来解决，正如我们在这项工作中所做的那样。此外，考虑到fMRI数据中头部运动体积的存在，图5中报告的高阶脑图似乎是稳健的(另见SI图14)。最后，我们强调，除了超复杂性指标外，我们的框架主要检测相干同步模式，而它大多忽略了退相干模式的影响，这在系统的整体动力学中是很重要的。未来的工作应该探索以更明确的方式处理数据中存在的退相干模式的替代方法。