在复杂系统的研究中, 组件间错综复杂的相互作用以及这些相互作用如何引发集体动力学行为一直是统计物理和复杂网络领域研究者的关注焦点^[1-6]. 通过分析系统动力学特性来推断其内部连接方式, 对于理解并控制复杂系统的整体功能至关重要. 因此, 利用网络节点的动力学行为来逆推网络结构成为了该研究领域的一个重要研究内容, 区别于仅考虑两个节点的传统低阶相互作用, 高阶相互作用则同时考虑了三个或更多节点之间所构成的局部集团共同对系统的作用^[7,8]. 研究表明, 高阶相互作用在生态系统^[9]、信息传播^[10]和疾病传播^[11]等研究中有重要意义. 目前, 基于非线性动力学理论的方法主要关注于从系统的集体动力学中研究节点间的直接相互作用. 这通常涉及首先根据实际数据构建合适的数学模型, 然后确定模型中的未知参数, 进而揭示网络的拓扑结构和节点间的非线性相互作用^[12-21]. 然而, 这些方法往往针对特定问题而设计, 只在特定条件下适用, 缺乏普适性. 为了克服这一局限性, 结合非线性网络动力学的无模型预测技术, 研究者们提出了利用外部驱动来实现网络连接重构的策略. 在这些无模型预测方法中, 通常采用统计关联计算(如皮尔逊相关性、斯皮尔曼相关系数、格兰杰因果关系等)来构建网络的功能连接^[12,22]. 然而, 值得注意的是, 皮尔逊相关系数仅适用于衡量符合正态分布数据的线性相关性, 对于非线性或非正态分布的数据, 它可能无法准确捕捉变量间的关系. 此外, 格兰杰因果性对数据中的异常值极为敏感, 处理长时间序列数据时容易出现偏差.

在揭示网络的连接结构时, 有效区分直接与间接连接至关重要. 例如, 可以通过分析系统受到的共同外部输入或网络中其他元素引起的去相关效应来研究间接相互作用^[23]. 虽然目前许多研究基于网络连边的统计关联提出了过滤间接连接的方法^[22,24], 但这些方法通常只在特定条件下适用(包括给定的网络尺寸、数据的稀疏性、数据的线性依赖性、噪声等条件), 而完全满足条件的案例在实际系统中很少出现. 基于非线性相位同步的概念^[23], 研究者提出了一种使用偏相关系数矩阵来减少间接连接影响的策略. 通过将偏相关系数矩阵与原始相关系数矩阵进行比较, 可以推断出网络的真实连接结构. 然而, 这种方法在节点数量较少的网络中效果较好, 并且要求系统处于近似同步状态, 这意味着它无法有效识别在完全非同步和完全同步状态下的相互作用. 因此, 尽管该方法在某些情境下是有益的, 其应用范围仍受限, 需要进一步的研究来扩展其适用性.

在探索细胞等较大尺寸网络的相互作用模式中, 文献[22]提出了一种全局静默方法. 首先, 研究者提出了全局响应矩阵G和局部响应矩阵S两种用于描述细胞组分间相互作用的矩阵. 全局响应矩阵$ G_{ij} = \dfrac{{\mathrm{d}}x_{i}(x_1, x_2, \cdots, x_{i-1}, x_{i+1}, \cdots)}{{\mathrm{d}}x_{j}} $捕捉了节点i对节点j变化的全局响应, 这种全局视角同时包含了直接和间接的两种相互作用, 无法进行两者的有效区分. 相比之下, $ S_{ij} = \dfrac{\partial{x_i}}{\partial{x_j}} $专注于局部响应, 通过数学变换“沉默”消除了间接效应, 而仅保留直接相互作用, 从而更准确地预测直接链接. 简而言之, G提供了一个包含所有相互作用的全面视图, 而S则通过消除间接效应而专注于直接相互作用的预测和理解. 然而值得注意的是, 文中消除间接效应时存在两个问题, 一定程度上限值了该方法的普适性. 1)在计算S时, 需要一个基于网络中扰动传播特性的近似条件, 即两个节点之间的响应主要由它们之间的最短路径决定, 从而忽略那些通过更长路径传播的间接效应, 而该近似条件所导致的误差必然将随着网络尺寸的减小而增大. 2)该间接作用消除法中需要对G进行求逆, 因此G可逆是这一方法的隐藏条件. 然而G矩阵的奇异性随着节点数的增加而显著增加. 因此, 该方法对于尺寸较小或较大的网络都不能很好地识别.

网络的拓扑结构往往可以通过节点之间的同步情况来判断. 当网络系统出现一般常见的同步情况时, 上述基于相关系数矩阵的大多数方法都是适用的. 当出现下列反常同步现象时, 上述基于相位同步的方法不再适用. 例如, 间接相连的节点间比直接相连的节点间提前出现的遥同步^[25,26], 耦合强度弱的一对节点比耦合强度强的一对节点优先出现的弱赢家相位同步^[27], 全同节点出现了同步与非同步混合的奇异态^[28,29]等. 文献[30]研究表明, 基于遥同步的思想, 利用一阶Kuramoto模型发现了网络中对称节点能够产生完全同步, 并在大脑功能网络结构的重构中得到了很好的验证. 需要指出的是, 复杂网络结构中的非对称节点也有可能产生同步, 因此该研究提出的同步是判断网络对称节点的充分非必要条件. 本文正是基于这一问题, 研究节点同步与网络结构识别的必要性条件.

近年来, 网络高阶相互作用的研究得到了快速发展, 比如同时考虑三个或更多节点所构成的局部集团的共同动力学作用^[31]. 由于传统的低阶相互作用受到其他高阶相互作用的影响, 从动力学的角度判断网络结构的存在性变得更为复杂, 因此网络高阶作用的出现进一步增大了揭示系统连接结构的难度. 本文将主要考虑以单纯复型为代表的高阶相互作用, 这对于预测系统行为、网络鲁棒性等具有重要意义. 本文提出了一种新方法, 用于推断耦合动力系统中的直接相互作用. 该方法基于同质化非线性动力学模型, 仅依赖于网络各节点的动力学特征, 能够精确地揭示节点间的直接相互作用, 而无需依赖统计依赖关系. 此外, 该方法还能有效识别网络中的对称节点. 我们还发现, 该方法同样适用于具有高阶相互作用的网络.

首先考虑单纯复形高阶相互作用下的Kuramoto相位节点的广义动力学模型及其广义高阶拉普拉斯矩阵^[32]. 具体而言, 随时间演化的动力学系统为

其中$ \theta_i \in [0, 2\pi) $是描述节点i的一维状态变量, 其动力学演化规则通常取$ f(\theta_i) = \omega_i $, $ g^{(D)}(\theta_i, \theta_{l_1}, \cdots, \theta_{l_D}) = \sin{ \Big(\displaystyle\sum\nolimits_{d = 1}^D \theta_{l_d} - d \theta_i \Big)} $; $ \lambda_i $表示系统第i阶的耦合强度; D表示系统所包含的最高阶数, d表示每一个阶数. 当$ D = 1 $时(一阶相互作用), 上述方程即为经典的Kuramoto模型, 当$ D > 1 $时, 系统不仅包括相互之间的低阶相互作用, 而且包含了更一般的高阶耦合作用, 由此动力学特性将会逐渐丰富. 图1展示了$ D\in[1, 3] $阶单纯复形高阶网络结构. $ a_{i, l_1, \cdots}^{(d)} $表示节点$ i, l_1, \cdots $之间的d阶邻接矩阵元, 其对应的d阶邻接矩阵$ {\boldsymbol{A}}^{(d)} $定义如下: 若存在一组d阶的相互作用$ \{i, l_1, l_2, \cdots, l_n, j\} $, 则$ a_{i, l_1, l_2, \cdots, l_n, j}^{(d)} = 1 $, 若不存在则$ a_{i, l_1, l_2, \cdots, l_n, j}^{(d)} = 0 $.

类似地, d阶相互作用的广义拉普拉斯矩阵可定义为

其中$ \delta_{ij} $为克洛内克函数,

当$ d = 1 $时, 则回归至一阶拉普拉斯矩阵$ L_{ij}^{(1)} $, 即:

首先分析一阶相互作用(传统低阶的)情况, 考虑包含相位延迟的Kuramoto节点:

其中$ \omega_i $为节点i的本征频率, $ \lambda_1 $为一阶相互作用的耦合强度, α为相位延迟项^[4]. 当节点趋近全同步时, 相位差较小, 因此上述方程在同步流形附近的线性化方程可表示为

其中, $ L_{ij}^{(1)} $即为一阶拉普拉斯矩阵((5)式).

在同步时, 不失一般性, 可以令$ \lambda_1 = 1 $以及$ \omega_i = 0 $. 此时, (7)式可以写为

考虑任意两个节点同步时, 可以得到关系式 $ \dot{\theta}_i = \dot{\theta}_j = \varOmega $, 即

其中$ \langle k \rangle ^{(1)} = N^{-1}\displaystyle\sum\nolimits_{j = 1}^{N}k_j^{(1)},\; \boldsymbol{1}=[1,1,\cdots,1]^{\rm T}$, 两个节点间的相位差则需满足下列方程:

我们发现, 对于Kuramoto-Sakaguchi (KS)节点所驱动的系统, 当系统中任意两个节点的相位与拉普拉斯矩阵之间满足关系式(10)时, 即发生同步. 上述关系式的矩阵形式为

其中$ {\left\langle {} \right. {\boldsymbol{k}} \left. {} \right\rangle } = {{{\mathrm{diag}}}}(k_1, k_2, \cdots, k_N) $是一个对角矩阵, 每个元素是对应节点的度; $ {\boldsymbol{\theta}}^{{\bf{\ast}}} $为

根据(11)式即可从节点相位时间序列中找到同步节点, 进而反推出网络的连接结构.

当仅考虑二阶而忽略一阶单形耦合时, KS动力学方程可表示为

此时的同步性由二阶拉普拉斯矩阵决定:

其中

即二阶拉普拉斯矩阵也可表示为

考虑任意两个节点同步时, 可以得到关系式 $ \dot{\theta}_i = \dot{\theta}_j = \varOmega $, 即

每个的相位可表示为

和一阶单形相互作用类似, 任意两个节点在同步条件下的相位差与拉普拉斯矩阵之间需满足如下关系:

由上述推导可知, 同时考虑一阶和二阶单纯复形耦合项时的动力学方程可以表示为

同样不失一般性, 设$ \omega = 0 $和$ \lambda_1 = \lambda_2 = 1 $, 则(20)式可以简化为

将拉普拉斯矩阵代入(21)式, 可得

为简便起见, 做如下定义:

考虑关系式$ \dot{\theta}_i = \dot{\theta}_j = \varOmega $, 依然有

于是, 每个节点的相位可表示为

于是得到与方程(19)类似的同步关系矩阵方程:

当系统包含从一阶到D阶单形的相互作用, 通过完全相同的方法可以得到关于广义拉普拉斯矩阵与系统各节点相位的方程, 同步条件为

定义:

则每个节点的相位可表示为

则对应的同步关系矩阵方程为

接下来证明Kuramoto-Sakaguchi系统的同步与网络结构的对称性有关. 一个由N个节点构成的网络, 其高阶拉普拉斯矩阵为$ {\boldsymbol{L}}^{(1{\text{-}}D)} $. 若能够找到一个双射$ \pi: {\boldsymbol{N}} \to {\boldsymbol{N}} $保持G的邻接关系, 则网络$ {\boldsymbol{G}}(N, {\boldsymbol{L}}^{(1{\text{-}}D)}) $有一个对称, 即G是一个自同构^[33]. 形式上, 这意味着存在一个置换矩阵$ {\boldsymbol{P}} = {\boldsymbol{P}}(\pi) $, 使得$ {\boldsymbol{P}}{\boldsymbol{L}}^{(1{\text{-}}D)}{\boldsymbol{P}}^{-1} = {\boldsymbol{L}}^{(1{\text{-}}D)} $. 如果自同构G的P与$ {\boldsymbol{L}}^{(1{\text{-}}D)} $是对易的, 那么$ {\boldsymbol{P}}{\boldsymbol{L}}^{(1{\text{-}}D)}{\boldsymbol{P}}^{-1} $操作保持原始网络邻接矩阵$ {\boldsymbol{L}}^{(1{\text{-}}D)} $的信息^[33]. 一般来说, 一个网络可以有m个自同构, 即这个网络有m组对称节点. 将系统在同步条件下的方程(9)两侧均左乘矩阵P, 可得

根据P和$ {\boldsymbol{L}}^{(1{\text{-}}D)} $的对易性质, 则有

与方程(28)对比可得到该系统存在交换对称性的方程为

(35)式表明: 通过左乘交换对称节点操作的交换矩阵P, 得到了交换后的拉普拉斯矩阵$ {\boldsymbol {PL}} $, 则一定有$ {\boldsymbol {PL}}{\boldsymbol{\theta}} = {\boldsymbol L}{\boldsymbol{\theta}} $. 又因为$ {\boldsymbol{PL}} $和L等价, 所以两者对应的θ也是相等的, 因此对称节点始终有相同相位. 另外, 由于上述变换未经过线性化近似, 故其判断的有效性将不依赖于相位延迟α的大小, 因此可以用任意给定的α来寻找任意系统的对称节点.

根据方程(35), 可以得出以下结论: 网络结构对称节点的相位在系统处于稳态后任意时刻均保持相同, 即任意的对称节点对i和j, 满足$ {\textit{π}} = {\textit{π}} (j\to i) $置换对称时, 存在结论$ |\theta_i-\theta_j| = 0 $. 由(35)式可以得出网络节点的同步与节点在结构上的对称是充分必要条件. 由此, 本文提出网络重构方法, 主要包括以下4个步骤(如图2所示).

1)待求系统网络节点KS方程同步动力学数值模拟.

对于由N个节点组成的系统, 由于不清楚原网络的单纯复形的阶数以及数量, 对该系统施加的KS方程如下:

式中取$ D = N-1 $, 即系统动力学包含了1到$ D = N-1 $阶的相互作用, D为最大可能阶次相互作用.

2)不失一般性, 采取频率全同节点并令所有节点的频率$ \omega = 0 $, 耦合强度$ \lambda = 1 $, 使用随机初值迭代方程.

3)根据全同节点的特点, 仅需取任意大小的耦合强度来克服随机初值的影响, 所有节点很快进入同步状态. 因此, 一定的暂态过后, 节点相位序列进入稳态, 可以获得每一个节点的相位时间序列.

4)根据(32)式计算拉普拉斯矩阵.

5)反推网络结构. 首先, 根据(36)式得到每个节点到达稳态后的的相位, 利用(32)式得到系统的广义拉普拉斯矩阵. 其次, 根据广义拉普拉斯矩阵的定义式(30), 就可以推导出一阶到D阶拉普拉斯矩阵从而得到每一阶的拉普拉斯矩阵, 即对应于每一阶单纯复形的相互作用.

接下来, 运用上述方法处理两类网络结构的问题: 1)以区分直接和间接连接为例, 识别传统的网络一阶相互作用; 2)识别网络一阶和二阶单纯复形混合的相互作用.

首先以直接和间接作用为例, 说明上述方法可以正确识别网络中一阶复形的相互作用. 以三个节点构成的系统为例, 在仅考虑一阶相互作用时, (31)式展开后可以得到其同步关系矩阵:

式中α为预先给定的参数值. 由于网络尺寸通常较大(即$ N > 3 $时), 只选取一个α得到的矩阵方程无法得到唯一的拉普拉斯矩阵解. 此时, 可以通过引入多个预设α而得到多组方程, 最终确定唯一的拉普拉斯矩阵解. 于是, (37)式的矩阵形式如下:

其中

(38)式左侧为元素已知的系数矩阵, 记为$ {{\boldsymbol{C}}_{{{n}}}} = {{\boldsymbol{\theta}}^{{\bf{\ast}}}{{(}}{{\alpha}}_{{n}}{{)}}} $. 记$ {{\boldsymbol{L}}^{{{(1)}}}} = \{L_{ij}|i = 1, 2\cdots, N; \ j >i\} $, $ {\boldsymbol{I}} $为单位列向量. 则(39)式可以简写为

当考虑1到D阶单纯复形相互作用时, (40)式可写为

在不同的$ \alpha_n $下进行数值模拟得到$ C_n $. 因此, 理论上可以通过寻找多个系数矩阵来求解拉普拉斯矩阵的所有矩阵元的数值.

例1　以三个节点组成的简单结构为例, 说明仅通过同步动力学求解网络结构的过程. 由方程(38)出发, 可以将网络拉普拉斯矩阵的各元素与每一对节点相位差的关系显式写为

图3展示了两种$ N = 3 $的网络结构在$ \alpha = 0.1 $时每一对节点间的相位差时间序列. 将图3(a)和图3(b)描述的系统稳定后任意时刻的节点相位(差)分别代入(42)式, 得到两个网络拉普拉斯矩阵为

从(43)式可知, 图3(a)和图3(b)描述了两种不同的网络结构. 因此, 这一方法有效地区分了直接和间接的作用.

例2　当网络尺寸增大时, 上述方法依然可以有效地识别网络中的对称节点, 并且对揭示网络中的遥同步现象具有重要的意义^[30], 结果如图4所示. 首先, 考虑一个由12个节点组成的网络, 其中节点1—节点10构成全连接网络, 节点11和节点12随机接入该网络, 如图4(a)所示. 图4(b)展示了节点11和节点12之间相位差的时间片段. 可以发现对于任意时刻, 相位差均满足$ |\theta_{11}-\theta_{12}|\equiv0 $. 对称节点具有相同的相位这一结论与该网络关 于节点11和节点12的交换对称群完全相符. 同样, 可以根据(38)式求出相应的拉普拉斯矩阵, 见附录A. 例2中, 实际上已经包含了一阶到D阶的耦合作用, 也就是可以进行$ {\boldsymbol L}^{(1)} $, $ {\boldsymbol L}^{(2)} $等的计算, 由于L的计算量太大, 本文忽略展示此时的L.

从网络结构的角度来看, 节点11和节点12之间没有直接连边, 而是通过中继全连网络间接相连. 此时它们之间产生的同步被称为遥同步现象, 这可以用更为简洁的结构图来说明. 比如, 图4(c)和图4(d)展示了网络结构的对称性对于遥同步现 象的研究具有重要意义. 以图4(c)所示的$ N = 5 $的星形网络为例, 在特定的参数选取情况下, 4个叶子节点(度$ k_{{\rm{leaf}}} = 1 $的节点)将先于其各自与中间节点(度$ k_{{\rm{hub}}} = 4 $的节点)发生同步. 我们发现正是因为该模型满足叶子节点的交换对称性, 即$ {\textit{π}}\; (1, i = 2:5) $, 这里的$ i \in [2, 5] $可以随意置换. 这与数值模拟的结果相符合(如图4(d)所示), 各叶子节点的相位差均为0, 而叶子节点与中心节点存在一个恒定且相同的相位差. 进一步通过上述方法所识别的网络拉普拉斯矩阵见附录A.

例3　本文在KS节点同步流形的附近, 采用线性近似的方法找到了同步与网络结构的关联性, 特别强调了相位延迟项α在网络识别中的重要作用. 以三个节点组成的链式网络为例, 说明线性近似的方法在较大范围内都是可行的, 结果如图5所示, 展示了各对节点相位差随着α的变化情况. 其中黑色虚线为根据方程(31)计算的理论解, 结果显示, 对称节点2和3之间相位差为$ \theta_{23} = 0 $(水平黑色虚线), 且与α的取值无关, 绿色的点线为数值模拟结果. 而对于非对称节点$ \theta_{12} = \theta_{13} = - {\alpha}/{3} $(倾斜黑色虚线是方程(31)的理论解), 分别与蓝色和橙黄色点线的数值模拟结果对应. 进一步可以发现对于任意的α, 对称节点的数值解和理论解均完全符合, 而对于非对称节点, 仅当α较小时数值解和理论解才较为相符.

本节展示上述方法可以有效区分不同阶次单纯复形的作用. 在图6(a)所示的网络结构中$ (N = 8) $, 同时存在一个一阶和二阶单纯复形混合作用. 在计算分析之前, 我们并不清楚网络单纯复形的阶数以及数量, 因此对该系统施加的动力学为

其中, $ D = 7 $为最高可能阶次. 根据网络对称节点的定义, 图6(a)所构成的群具有交换对称的关系

即图6(a) 中仅节点1和节点2为对称节点. 为了简单起见, 将节点1作为参考节点, 计算其他节点与节点1之间的相位差, 如图6(b)所示. 其中仅 有$ |\theta_2 - \theta_1 | $在系统达到稳定之后恒定为0, 而其他所有节点与节点1之间的相位差均收敛为非零 值. 因此得到了该系统中仅节点1和节点2对称的结论.

通过(28)式, 可以得到预测的广义拉普拉斯矩阵L:

接下来, 根据定义式(30)对上式L进行分解, 由此得到一阶到D阶拉普拉斯矩阵如下:

于是, $ {\boldsymbol{L}}^{(1)} $即重构的网络, 而$ {\boldsymbol{L}}^{(2)} $为可能存在的二阶单纯形作用, 这与图6是完全匹配的.

本文提出了一种基于同步动力学来重构网络单纯复形结构的方法. 通过分析发现KS模型的时间序列数据与网络拓扑结构之间存在精确的解析联系. 这一重要发现使得我们能够准确识别并重现网络的结构. 该方法不仅具有较强的普适性, 而且对于包含任意阶单纯复形相互作用的网络结构(即d阶单纯复形), 其识别效能同样显著. 在规模为$ N = 2 $到$ N = 10 $的网络中, 成功地辨识了网络的结构. 此外, 还探讨了该解析关系在识别网络对称节点方面的应用价值, 证实了网络结构的交换对称性在动力学上表现为完全同步状态, 并通过数值模拟验证了这一理论. 在典型的星形网络模型上进行了实证, 进一步的研究揭示了遥同步现象与网络结构对称性之间的依赖关系. 这样的研究为理解网络动力学的同步机制提供了新的理论视角. 基于多维时间序列进行高阶网络的识别, 目前还处于起步阶段, 本文的方法和其他方法的对比分析非常值得进一步的研究. 今后, 更应该关注更大尺寸网络上高阶相互作用的识别, 尤其是实际复杂网络上的单纯复形或者有向复形相互作用的动力学模型构建和实证研究.

图4中, 节点11和节点12通过全连网络间接相连网络结构所识别的(一阶和二阶单纯复形)拉普拉斯矩阵如下:

星形网络(图4(c))所识别的拉普拉斯矩阵如下:

0表示所有元素均为0的矩阵.