系统弹性是复杂系统的重要特征之一，描述系统受到干扰、发生错误和环境变化时，调节自身以保持基本功能的能力^[1]。人们往往是在系统失去弹性、引发毁灭性后果时才开始采取措施。但此时已造成的损失不可挽回，并且系统恢复弹性非常缓慢。因此，在系统失去弹性、产生严重损失之前，准确度量系统状态和测量系统弹性，从而及时采取措施提高系统弹性是防止系统崩坍的重要方法之一。

系统崩塌对生态系统稳定^[2]和社会经济发展^[3]造成广泛影响，但导致系统丧失弹性的事件往往难以预测。其根源为多数系统弹性的研究都关注由少数成分构成的低维模型，而不是更符合实际的、由多属性构成的多维复杂模型。2016年，文献[4]提出复杂系统弹性预测模型，通过将多维复杂系统映射至一维空间，有效预测复杂系统弹性，弥补了多维复杂系统弹性和一维动力学系统弹性间的理论差距。该模型提出后，受到广泛专注。文献[5]发现文献[4]对模型所做限制过于苛刻，并拓展了模型的应用条件。文献[6]基于弹性预测模型，介绍了复杂网络“质量”中心的概念，证明系统的相变和状态演化以“质量”中心为活动特征。文献[7]研究了复杂系统弹性的抗噪能力，发现噪声可以妨碍或提高临界跃迁。

在网络环境中，每个节点的状态不仅受其一阶邻居状态的影响，还受其高阶邻居状态的影响。高阶邻居信息被运用于多个领域。在节点相似性领域，基于两个节点最近邻相似则这两个节点相似的概念，文献[8]提出了一种网络中顶点对结构相似性的度量方法。针对推荐系统中的协同过滤算法，文献[9-10]通过考虑二阶相关性降低主流偏好的影响，分别提出了改进后的协同过滤(MCF)算法和有向二阶协同过滤(HDCF)算法，有效提高推荐算法的准确性和多样性。然而复杂系统弹性预测模型^[4]中的算子仅考虑节点最近邻信息，并没有考虑节点二阶邻居的影响。同时，网络中由于节点和边之间存在不同的联结方式，就具有不同的拓扑结构和特性、行为，因而需要用不同的数学模型来描述^[11]。具有代表性的网络模型包括规则网络、随机网络、小世界网络和无标度网络^[12]。这4类网络模型作为复杂网络理论的基础，被用于研究Internet网^[13-14]、病毒传播^[15]和电力交通^[16-17]等系统。

受高阶信息应用的启发，本文提出基于节点二阶邻居信息的弹性预测模型，并在平均度分别为2、4、6、8、10的BA无标度网络和WS小世界网络上验证改进后模型的有效性，进而探索网络拓扑结构对改进后模型效果的影响。实验首先对每个人工网络重复进行不同程度的随机删节点、删边、改变权重这3种干扰。然后运用基因调控系统动力学方程，计算系统稳定时，干扰后网络中各节点的状态，并用原始模型和基于节点二阶邻居信息的复杂系统弹性测量模型分别预测系统弹性。最后，为量化改进后模型的测量复杂系统弹性效果，本文定义了指标系统弹性测量精度$W$。实验结果表明，基于节点二阶信息的弹性预测模型在不同平均度的两类网络中均有效。和原始模型相比，BA无标度网络和WS小世界网络的系统弹性测量精度分别最多提高79.89%和59.53%。对于同规模同平均度的网络，改进后模型在BA无标度网络中的效果优于WS小世界网络。整体而言，网络平均度越小，改进后模型的效果越明显。本文的工作进一步证明考虑二阶信息是一种有效提高系统弹性测量准确性的方法。

1.   基础理论

1.1   弹性预测模型

由$N$个组件构成的系统，其中所有节点都具有活性特征，即${\boldsymbol{x}} = {({x_1},{x_1}, \cdots ,{x_N})^{\rm{{\rm T}}}}$。该系统的动力学方程依赖于以下耦合非线性方程^[18]：

式中，$F({x_i})$描述每个组件自身的动力学特征；$\sum\limits_{j = 1}^N {{a_{ij}}G({x_i},{x_j})} $描述节点$i$和它邻居节点的相互作用；${a_{ij}}$表示邻居节点$j$对目标节点$i$的影响；函数$G({x_i},{x_j})$意味着节点$i$和节点$j$之间的动力学特征。当选取合适的$F({x_i})$和$G({x_i},{x_j})$时，式(1)可用于仿真不同的系统，如单元系统^[19]、生态系统^[20-21]和社会系统^[22]。

为了将多维方程映射为一维方程，并有效预测系统弹性，文献[4]提出一种算子，形式如下：

式中，单位向量${\boldsymbol{l}} = {(1,1, \cdots ,1)^{\rm{{\rm T}}}}$；${\boldsymbol{y}} = {({y_1},{y_2}, \cdots ,{y_N})^{\rm{{\rm T}}}}$，${y_i}$表示和节点$i$相关的标量，例如节点$i$的活性${x_i}$；${\boldsymbol{A}}$为加权邻接矩阵。算子$L$平均目标节点的所有邻居节点标量值作为输出，${\boldsymbol{y}} = {({y_1},{y_2}, \cdots ,{y_N})^{\rm{{\rm T}}}}$作为输入。对于计算过程，式(2)可以写成：

式中，$s_j^{{\rm{out}}} = \sum\limits_{i = 1}^N {{a_{ij}}} $。如果${y_j}({x_i}) = G({x_i},{x_j})$，当度不相关时，那么节点$j$的均值就独立于节点$i$。换言之，假设$i$的最近邻均值和其他所有节点的最近邻均值相同，则式(1)中由$i$的所有邻居$j$作用之和构成的交互项可以写成：

式中，$s_i^{{\rm{in}}}{\rm{ = }}\sum\limits_{j = 1}^N {{a_{ij}}} $；${\langle {y_j}({x_i})\rangle _{nn}}$表示节点$i$的邻居节点状态${y_j}({x_i})$的均值。此外，${\boldsymbol{x}} = {({x_1},{x_2}, \cdots ,{x_N})^{\rm{{\rm T}}}}$，$L(G({x_i},{\boldsymbol{x}}))$表示所有最近邻节点$j$的均值。将式(4)代入式(1)，则：

对式(1)基于平均场理论，有$L(G({x_i},{\boldsymbol{x}})) \approx $ $G({x_i},L({\boldsymbol{x}}))$。式(5)可以写成：

向量的形式为：

式中，$\langle \circ \rangle $表示Hadamard卷积^[23]。将式(2)应用于式(7)的两侧，则有：

式中，基于前文提到的平均场假设，有$L(F(x)) \approx $ $F(L(x))$和$L({{\boldsymbol{s}}^{{\rm{in}}}} \circ G(x,x)) \approx L({{\boldsymbol{s}}^{{\rm{in}}}})G(L(x),L(x))$。

邻居活性${x_{{\rm{eff}}}}$可以表示为：

并且，最近邻加权度${\beta _{{\rm{eff}}}}$表示为：

然后，式(8)可被写为：

它将式(1)表示的多维系统映射为有效的一维方程，该方程以一个控制参数${\beta _{{\rm{eff}}}}$为特征。

1.2   网络构造算法

不同拓扑结构的网络具有不同的性质，本节对WS小世界网络和BA无标度网络这两种网络的构造算法进行简单介绍：

1) WS小世界模型构造算法$G(N,K,p)$^[24]：

① 给定一个含有$N$个点的环状最近邻耦合网络，其中每个节点都与它左右相邻的各$K/2$个节点相连，$K$是偶数。

② 以概率$p$随机的重新连接网络中原有的每条边，即每条边的一个端点保持不变，另一个端点取网络中随机选择的一个节点。其中规定不得有重边和自环。

2) BA无标度网络模型构造算法$G(N,m)$^[17]：

① 增长：从一个具有${m_0}$个节点的连通网络开始，每次引入一个新的节点并且连到$m$个已存在的节点上，这里$m \leqslant {m_0}$。

② 优先连接：一个新节点与一个已经存在的节点$i$相连接的概率${\mathit{Π} _i}$与节点$i$的度${k_i}$之间满足如下关系：

对于不同类型网络特征的描述可参考文献[16]。

2.   模型及实验设计

2.1   实验数据及模型构建

实验数据集包括两种不同类型的网络：BA无标度网络和WS小世界网络，每种类型生成平均度为2、4、6、8、10的5个网络，网络规模均为100，人工生成10个网络作为实验数据集。

考虑节点的二阶邻居信息，本文改进了原模型的算子。基于节点二阶信息的算子，形式如下：

式中，${{\boldsymbol{A}}^2}$描述目标节点和它的二阶邻居节点之间的相互影响。另外，$\boldsymbol{y} = {({y_1},{y_2}, \cdots ,{y_N})^{\rm{{\rm T}}}}$中的${y_i}$是一个有关目标节点$i$的标量。

则加权平均邻居活性${x_e}$可以表示为：

加权平均度${\beta _e}$表示为

式中，${\boldsymbol{x}} = {({x_1},{x_2}, \cdots ,{x_N})^{\rm{{\rm T}}}}$为节点活性的向量；$s_i^{{\rm{in}}} = \sum\limits_{j = 1}^N {{a_{ij}}} $且${{\boldsymbol{s}}^{{\rm{in}}}} = (s_1^{{\rm{in}}},s_2^{{\rm{in}}}, \cdots ,s_N^{{\rm{in}}})$。

利用新的算子形式$Z({\boldsymbol{y}})$替代原始算子$L({\boldsymbol{y}})$，多维系统动力学通式被映射为如下的一维方程：

上式等于0时，解析解表示系统可能的状态值${x_a}({\beta _e})$。且当式(16)对变量${x_e}$求导小于0时，保证线性稳定性，这时可以得到系统稳态的数值分析结果${x_a}$。

2.2   实验设计

本文计算每个节点状态时，基于基因调控系统的动力学方程Michaelis-Menten为^[18]：

式中，右侧第一项描述每个细胞的自动力学；第二项为细胞间活动，其中希尔指数$h$代表基因调控中的合作水平^[18]。本文设定$B = 1$，$q = 1$，$h = 2$，为大肠杆菌转录网络的动力学方程，对两种类型的10个网络分别进行干扰实验，探索基因调控系统的弹性变化。发现当细胞承受足够强的干扰时，会从弹性状态转为死亡态。具体实验过程如下：首先，对网络进行3种不同方式不同程度的扰动，包括随机删除${f_n}$比例的节点、随机删除${f_l}$比例的边、改变${f_w}$比例的全局权重。这3种扰动分别模拟基因剔除、抑制相互调节作用和环境变化。初始所有节点的状态设为${x_0} = 2$，并依赖于动力学方程计算系统收敛时的节点状态。对于每一种的不同程度扰动，执行100次的移除。结果表明，细胞可以承受一定程度的干扰，而当扰动足够强时细胞从弹性状态($x > 0$)转化为死亡态($x{\rm{ = }}0$)。对每一次干扰，无论何种扰动类型，过渡点(网络的平均状态$\langle x\rangle $首次变为零)的值都不同。然后，基于上文提到的映射过程，式(17)重写成式(16)一维的形式为：

此过程将多维方程公式映射到$\beta $-空间，每个网络受到一定程度的某种干扰时，生成一个新网络。同一个网络依次进行不同方式不同程度的干扰，会生成$n$个网络，将$n$个网络映射至$\beta $-空间，每个网络映射为一个点。最后，本文定义系统弹性测量精度，量化基于节点二阶信息的弹性测量模型的效果。

2.3   衡量指标

本节定义指标系统弹性测量精度$W$来对比提出的改进后新模型和原始模型，即考虑最近邻信息模型的效果。$W$可表示为：

式中，$\sum\limits_{l = 1}^n {{{({x_{{\rm{al}}}} - {x_{{\rm{el}}}})}^2}} $表示式(18)得到的系统状态的数值分析结果${x_{{\rm{al}}}}$和式(14)获得的实证数据集映射结果${x_{{\rm{el}}}}$之间的总误差。$\sum\limits_{l = 1}^n {{{({x_{{\rm{al}}}}^\prime - {x_{{\rm{el}}}})}^2}} $是原始模型的总误差。当$W$大于0时，说明考虑节点二阶邻居信息的模型可比原始模型更精确地测量系统弹性。

3.   实验结果

基于上文所提出的改进后模型，依赖基因调控系统动力学方程对同规模不同类型网络的弹性测量效果进行对比。对于规模为100平均度为2的两种类型网络，实验结果如图 1所示，其中横、纵坐标分别为入度均值${\beta _e}$和状态均值${x_e}$，深色点为原始模型的映射结果，浅色点为基于节点二阶信息弹性预测模型的映射结果。实线为由式(18)推导得出的基因调控系统弹性模式。图 1a为网络平均度为2的BA无标度网络的映射图。结果表明，改进后模型的映射结果(浅色点)，相比原始模型的映射效果(深色点)，更贴近弹性模式(实线)。图 1b为网络平均度为2的WS小世界网络的映射结果，也可得出相同的结果。

图  1  网络的弹性映射结果对比

下载: 全尺寸图片 幻灯片

为进一步量化新模型效果和探索不同网络构型对改进后模型效果的影响，基于上文提出的系统弹性测量精度$W$，计算规模为100的不同平均度的BA无标度网络和WS小世界网络的弹性测量精度，如图 2所示。

图  2  不同平均度不同类型网络的系统弹性测量精度

下载: 全尺寸图片 幻灯片

由图 2可知，运用原始模型和改进后模型，通过对原始网络进行多次不同程度不同类型的干扰，计算系统弹性的理论解析值和实验映射值的平均误差，得到系统弹性测量精度。例如，针对$\langle k\rangle $为2的BA无标度网络，每删一次节点生成一个新的网络，不同比例删节点会产生多个网络，删边和改权重类似。3种干扰方式，每种干扰方式仿真100次，因此会产生多个网络，再对多个网络进行映射，计算映射后的弹性测量精度，视为$\langle k\rangle $为2的BA无标度网络的弹性测量精度。如图 2所示，实验结果表明同规模同类网络中，网络平均度越小，基于节点二阶信息的模型越适用。同时，针对同类型平均度相同的网络，改进后模型在BA无标度网络上的效果优于WS小世界网络。其中，两类网络平均度为2时，系统弹性测量精度$W$提高最明显，BA无标度网络为79.89%，WS小世界网络为59.53%。

4.   结束语

准确测量系统弹性是维持复杂系统稳定性的关键，2016年文献[4]提出通过将多维复杂系统映射至有效一维动力学，准确预测系统弹性的模型。基于节点二阶信息，本文对模型进行改进，并在同规模不同平均度的BA无标度网络和WS小世界网络上进行实验，验证了改进后模型的有效性。同时，探索了网络拓扑结构对改进后模型效果的影响。基于本文提出的系统弹性测量精度指标$W$，度量新模型的改进效果。结果表明对于10个网络，改进后模型均可更精确地测量系统弹性。其中网络平均度为2时，系统弹性测量精度$W$提高最明显，BA无标度网络为79.89%，WS小世界网络为59.53%。在同规模同类网络中，网络平均度越小，基于节点二阶信息的模型越适用。同时，针对同类型平均度相同的网络，改进后模型在BA无标度网络上的效果优于WS小世界网络。

本文通过考虑节点二阶信息，改进了复杂系统弹性测量模型，有效提高了系统弹性测量精度，为有效度量复杂系统弹性状态提供支持，也为设计弹性系统、防止系统崩塌提供科学依据。与此同时，本文研究也存在诸多不足，如不同系统依赖于不同的动力学方程，在其他系统中，网络特征对系统弹性度量是否有相同的影响。同时，不同特征的网络弹性预测结果差异的原因也有待进一步研究。