集体动力学是理解自然界中复杂耦合系统自组织行为的核心范式. 这一范式的普适性源于一个基本事实: 有相互作用单元组成的群体, 常常涌现出同步^[1-3]、奇异态^[4-8]、波传播^[9,10]、斑图形成^[11,12]等集体动力学行为. 这一观点得到了横跨物理、生物和工程领域的复杂系统数十年理论与实验研究的有力支持: 从萤火虫的节律性闪光^[13,14]、心脏起搏细胞的集体放电^[15,16]、到电网的协同运行^[17,18]以及大脑中的大尺度振荡活动^[19,20], 无不体现了一个普适原理, 即局部动力学单元通过自组织协同, 能够涌现出维持系统功能的宏观有序态. 理解此类涌现动力学如何产生、维持并最终退化, 不仅具有基础科学意义, 也在神经形态计算、神经系统疾病诊断等应用领域具有重要价值.

然而, 现实中的复杂系统往往并非由完全健康的单元构成. 在生物衰老、器件老化的影响下, 部分单元可能逐渐丧失其内在动力学活性, 从自激振荡态退化为静息稳态, 从而转变为“非活性”组分. 比如, 在阿尔茨海默病和帕金森病等神经退行性疾病中, 大量神经元失去自发放电能力, 进入静息态, 不再参与大脑网络的节律性协同^[20]; 在心脏系统中, 衰老过程伴随起搏细胞数量减少及关键离子通道表达下调, 导致窦房结整体起搏活性显著降低^[21]. 当失活动力学单元的比例超过某一临界阈值时, 整个网络可能突然丧失宏观振荡能力, 从集体活跃态塌缩至全局静息态. 这一现象被称为老化转变(aging transition), 最早由 Daido 和Nakanishi^[22]提出, 他们在全耦合网络中研究了自激与非自激节点共存所引发的集体动力学塌缩. 振荡死亡^[23-28]由耦合诱导的稳定性改变所致, 而老化转变是一种由组分异质性引发的集体失活, 指因部分单元丧失内在动力学活性而引发的集体振荡塌缩现象, 与生物老化过程中的动力学无直接关联. 此后十余年, 老化转变的研究在理论与机制层面不断深化^[29]. 早期工作聚焦于全局耦合系统, 系统刻画了自激与非自激单元混合比例如何决定集体振荡的存续^[22,30-32]. 随后, 关于老化转变的研究被扩展至更贴近现实的复杂拓扑结构. 研究表明, 加权耦合网络被证明可通过调节连接强度分布显著调控临界失活比例, 并导出其与平均权重的解析关系^[33]; 而在具有距离依赖耦合的空间嵌入神经元网络中, 老化转变甚至呈现出爆炸性特征, 表现为序参量的不连续塌缩与滞回曲线, 揭示了一级相变类动力学行为^[34]; Tanaka等^[35]发现, 在无标度Kuramoto振荡器网络中, 老化转变的临界点强烈依赖于失活节点的选择方式: 当节点随机失活时, 低度节点因其数量优势, 对维持集体同步具有决定性作用; 而针对低度节点的选择性失活则会显著降低网络的动力学鲁棒性. 这一反直觉现象表明, 网络的动力学鲁棒性并非仅由高度枢纽决定, 而是在特定条件下依赖于大量弱耦合单元的集体协同. 大量研究表明, 老化转变及其鲁棒性具有广泛的普适性. 这主要体现在3个方面: 其一, 多种动力学模型均可呈现此现象, 包括经典的Stuart-Landau振荡器^[36-43]、量子相干体系^[44,45]及各类神经元模型^[46-49]等; 其二, 其行为在不同网络拓扑中均被观察到, 如全局耦合^[46,47]、小世界^[41,50]与无标度网络^[36,49]等; 其三, 多种耦合机制^[39,42,43], 如扩散耦合^[36,45]、时滞耦合^[37,51]以及神经元的电/化学突触耦合^[50]等, 均对其产生关键调控作用.

然而, 传统对老化转变的研究大多建立在“节点间仅存在成对相互作用”这一简化假设之上. 尽管这构成了复杂网络理论的基石, 但在诸多真实系统中, 尤其是大脑这样的高阶生物网络, 单元间的协同往往不局限于两两之间, 而是涉及3个或更多节点之间的群体效应, 即高阶相互作用^[52-55]. 例如, 神经元之间不仅通过突触形成双向连接, 还可能通过电突触、星形胶质细胞介导的调节、多突触回路等多种机制形成三体乃至更高阶的协同^[56,57]. 这种高阶相互作用已被证明能显著改变系统的同步能力^[58]、诱导出丰富的时空动力学斑图^[8,59]. 因此, 我们有理由推断, 高阶相互作用也可能对老化转变产生深刻影响.

目前, 关于高阶相互作用如何影响老化转变的研究刚刚兴起, 相关文献尚在积累之中. 目前仅有少数工作对此进行了初步探索. 近期, Sharma等^[60]在加权无标度网络和二维单纯复形上研究了Morris-Lecar神经元的老化转变, 发现Morris-Lecar神经元在无标度网络中的老化转变表现为一阶相变, 高阶相互作用会降低网络的鲁棒性, 并消除中等耦合强度下的部分振荡态. 几乎同期, Biswas 和Seth^[61]则基于Rulkov映射研究了全局耦合网络中高阶相互作用对老化转变的调控, 发现高阶耦合可平滑或延迟爆炸性转变, 其效果依赖于相互作用的阶数和强度. 这些发现共同表明, 高阶相互作用的动力学效应高度依赖于具体的神经元模型与网络结构.

基于上述背景, 本研究旨在从网络动力学的角度, 探究生物神经系统在局部功能损伤下的集体响应与恢复机制. 为此, 我们在全局耦合的FitzHugh-Nagumo (FHN)神经元网络中, 同时引入描述直接连接的成对相互作用(1-单纯形)和刻画群体协同的高阶相互作用(2-单纯形), 而且这里的2-单纯形采用的是一种非线性耦合. 通过逐步改变网络中可兴奋神经元(即功能“损伤”单元)的比例, 并正向(模拟损伤累积)与反向(模拟损伤修复)扫描系统状态, 我们深入揭示了系统的“损伤耐受性”与“恢复滞后性”. 研究发现: 仅当存在成对相互作用时, 系统表现出双爆炸性老化失活与清晰的滞回现象, 表明系统从振荡状态到死亡态的损伤过程与反向恢复过程并不可逆, 存在显著的恢复滞后——即系统“康复”所需的条件比“受伤”更为苛刻. 高阶相互作用的引入则会改变这一特性, 诱发单爆炸性转变, 并在与成对耦合协同作用时, 有效抑制滞回行为, 表现为滞回环随高阶耦合强度增强而收窄乃至消失, 同时可诱导出丰富的过渡态. 这些结果表明, 高阶相互作用不仅拓展了系统对损伤的耐受边界, 也调控着其恢复路径的滞回特性, 从而在维持网络动态鲁棒性中扮演了与成对相互作用不同却至关重要的角色. 本研究为理解复杂生物网络中功能损伤与恢复的非平衡动力学提供了新的机制见解和理论依据.

为探究高阶相互作用对神经元网络老化转变的影响, 本文构建了一个由$ N=300 $个FHN神经元组成的全局耦合高阶网络. 模型的核心在于同时考虑了神经元之间的成对(两体)相互作用与高阶(三体)相互作用.

在本模型中, 网络拓扑由两类连接结构定义. 模型中的成对相互作用(1-单纯形)基于全局耦合假设, 即任意两个神经元之间均存在连接. 高阶相互作用(2-单纯形): 若3个神经元两两之间均存在成对连接, 则它们构成一个2-单纯形, 代表三体相互作用. 本文不考虑比2-单纯形更高阶的相互作用. 为直观展示, 图1给出了一个由10个节点构成的高阶网络示意图. 图中每条边代表成对相互作用. 图中红色高亮部分展示了其中一个2-单纯形(三体相互作用), 为简化图示仅展示了一个2-单纯形, 它包含3个节点i, j, k. 在该单纯形中, 本文考虑的高阶相互作用为: 节点j和 k 联合驱动节点i, 对应方程(1)中的$ (x_{j}^{2}{x}_{k}-x_{i}^{3}) $项作用于节点 i 的动力学方程. 需要说明的是, 一个完整的2-单纯形包含6个可能的高阶驱动方向(每个节点可被其余两个节点联合驱动).

耦合系统的动力学由以下方程描述:

其中, $ {x}_{i} $和$ {y}_{i} $$ \left(i=1, 2, \cdots , N\right) $分别为第$ i $个FHN神经元的膜电位与恢复变量; $ \varepsilon $为时间尺度分离参数$ \left(\varepsilon \ll 1\right) $, 确保$ {x}_{i} $与$ {y}_{i} $演化速率不同, 本研究中, $ \varepsilon $取值为0.005; $ {k}_{\text{p}} $和$ {k}_{\text{h}} $分别为成对相互作用与高阶相互作用的耦合强度; $ \langle {k}^{1}\rangle $和$ \langle {k}^{2}\rangle $分别为网络的一阶平均度(平均成对连接数)与二阶平均度(平均参与的三体相互作用数), 用于对耦合项进行归一化. 由于本文采用全局耦合网络, 任意两个神经元之间均存在成对连接, 因此一阶平均度即为网络中每个节点的邻居数量, 也就是说$ \langle {k}^{1}\rangle =N-1 $. 对于高阶相互作用, 任意3个神经元两两相连即构成一个2-单纯形. 在全局耦合网络中, 每个节点参与的三体相互作用数量为从其余$ N-1 $个节点中有序选取两个的组合数, 因此二阶平均度为$ \langle {k}^{2}\rangle = (N-1)(N-2) $. 在没有任何耦合, 即$ {k}_{\text{p}}={k}_{\text{h}}=0 $时, 这些相互独立的FHN 神经元的动力学由参数$ {a}_{i} $决定. 当$ \left| {a}_{i}\right| \lt 1 $时, 神经元处于自振荡状态; 当$ \left| {a}_{i}\right| \gt 1 $时, 神经元处于可兴奋态, 即需要外部刺激才能产生神经脉冲. 在本文中, 对于可兴奋神经元, 取$ {a}_{i}=1.3 $; 对于自振荡神经元, 取$ {a}_{i}=0.7 $. 需要特别指出的是, 本模型在高阶相互作用项中引入了非线性耦合机制$ (x_{j}^{2}{x}_{k}-x_{i}^{3}) $^[62,63], 这与以往研究中常采用的线性耦合形式有本质不同. 线性高阶耦合在数学上可重新表述为等效的加权成对耦合, 而本文采用的非线性形式则能体现三体相互作用的不可约性, 是研究高阶动力学效应的更合适方式.

通过调节系统中初始可兴奋神经元占神经元总数的比例$ {P}_{\text{a}} $, 控制网络中可兴奋神经元与自振荡神经元的数目, 以研究系统在损伤梯度下的老化转变现象. 为清晰说明这一调控过程, 图2展示了一个示例. 取$ N=50 $个神经元为例, 其中图2(a)表示每个神经元的参数$ {a}_{i} $随时间的变化: 红色对应$ {a}_{i}=0.7 $(自振荡神经元), 蓝色对应$ {a}_{i}=1.3 $(可兴奋神经元). 图2(b)则表示可兴奋神经元比例$ {P}_{\text{a}} $随时间的变化. 初始时$ {P}_{\text{a}}=0 $系统中所有神经元均处于自振荡状态. 随后, 每经过一个固定时间间隔$ T $, 随机选择系统中5个自振荡神经元将其变为可兴奋状态, 直至$ {P}_{\text{a}}=1 $, 即所有神经元均处于可兴奋态. 为探究系统可能存在的双稳态及滞回行为, 进一步进行反向扫描: 每经过一个时间间隔$ T $, 将正向扫描过程中在此间隔内被失活的5个神经元恢复为自振荡状态, 直至系统中所有神经元恢复为自振荡态. 正向与反向过程保持完全对称, 以确保扫描路径的可逆性, 从而有效揭示系统在“损伤累积”与“损伤修复”过程中的滞后响应.

采用二阶龙格-库塔法对动力学方程组(1)和(2)进行数值积分, 时间步长设为$ \text{d}t=0.002 $. 本文中, $ N=300 $. $ T=300 $, 对应$ 1.5\times {10}^{5} $个时间步. 在每个$ T $的整数倍的时刻有4个神经元被“失活”或“激活”. 每个$ T $开始时的初始状态为上一个$ T $的末状态. 为消除瞬态影响, 舍弃每个$ T $时间段内前$ 1.0\times {10}^{5} $步的数据, 并在后续每间隔50个步长记录一次系统状态, 用于分析系统的稳态行为. 为确保系统演化满足准静态条件, 时间间隔$ T $的取值需保证系统在每次状态切换后充分驰豫至稳态, 每次切换的神经元数量$ M $需足够小以避免突变扰动. 经数值验证, 在$ T=250 $, 300, 350, 400及$ M=1 $—6的范围内, 本文的主要结论(爆炸性转变、滞回行为等)均保持稳定, 不依赖于$ T $和$ M $的取值. 对正向或反向扫描的每个$ {P}_{\text{a}} $(对应于一个特定的时间周期$ T $), 计算各个神经元的时间序列的振幅:

式中$ {x}_{i, \max } $和$ {x}_{i, \min } $分别表示变量$ {x}_{i} $的最高值与最低值. 如果$ {x}_{i, \text{amplitude}} $超过某个阈值, 这个神经元就处于振荡态, 否则就处于静息态. 在本文模拟中, 这个阈值取为2.5. 由于FHN神经元在振荡态时膜电位峰峰值与静息态时波动幅度存在显著量级差异, 因此该阈值在较宽范围内(如1.5—3.0)的选取对振荡神经元比例的统计结果无实质性影响. 系统状态的量化指标为振荡态神经元数量$ {N}_{\text{A}} $占神经元总数的比例:

当所有神经元都处于振荡态时$ {P}_{\text{N}}=1 $, 在完全失活(神经元全部振荡死亡)时$ {P}_{\text{N}}=0 $, 而部分失活(也是部分振荡)时, 则有$ 0 \lt {P}_{\text{N}} \lt 1 $.

作为基准, 当系统中既无成对耦合也无高阶相互作用时(耦合强度均为0), 系统表现为连续且可逆的转变(图3(a)). 随着可兴奋神经元比例$ {P}_{\text{a}} $从0增至1, 振荡神经元比例$ {P}_{\text{N}} $从1平滑下降至0; 反向扫描时, 系统沿相同路径返回, 未出现滞回现象. 这表明在无耦合作用下, 系统的动力学状态完全由各神经元的本征属性线性叠加决定. 值得注意的是, 图3曲线为11次独立随机重复实验的典型结果. 由于全局耦合网络的对称性及准静态变化过程, 不同随机种子下的曲线高度重合, 误差棒极小.

保持高阶耦合强度为零, 仅引入成对相互作用后, 系统的转变行为发生根本性变化, 呈现出爆炸性特征与滞回现象. 随着$ {P}_{\text{a}} $的增大或减小, $ {P}_{\text{N}} $在某一临界点发生突然的、不连续的跃迁(图3(c), (d), (g)), 即发生了爆炸性转变. 而且正向与反向扫描的转变临界点并不重合, 两条曲线分离形成滞回环, 表明系统状态具有路径依赖性, 即“损伤累积”与“损伤修复”并非可逆过程. 对比图3(c), (d), (g)可以发现, 随着成对相互作用强度$ {k}_{\text{p}} $的增大, 爆炸性转变从“不彻底”趋于“彻底”. 这里的“彻底”表示在完全振荡态和完全静息态之间的突然转变. 在$ {k}_{\text{p}} $较小时(图3(c), $ {k}_{\text{p}}=0.35 $), 系统从完全振荡态跃迁至一个部分振荡态($ 0 \lt {P}_{\text{N}} \lt 1 $), 随后$ {P}_{\text{N}} $才平滑降至0; 而当$ {k}_{\text{p}} $足够大时(图3(d), (g)), 系统直接从完全振荡态爆炸性地跃迁至完全振荡死亡态. 同时, $ {k}_{\text{p}} $增大时爆炸性转变的临界点显著左移. 这表明更强的扩散耦合非但未能提升网络的鲁棒性, 反而降低了系统发生老化转变的损伤阈值. 其物理机制在于强耦合加剧了失活动力学的传播效率, 使得局部失效更容易引发全局崩溃.

当系统中仅存在高阶相互作用时($ {k}_{\text{h}} \gt 0, {k}_{\text{p}}=0 $), 系统表现出单爆炸性转变(图3(b)). 在正向扫描中, $ {P}_{\text{N}} $在$ {P}_{\text{a}} $达到临界值时发生爆炸性下跌; 然而, 在反向扫描中, $ {P}_{\text{N}} $从0到1的变化却是连续平滑的, 未出现突然跃迁. 这表明高阶相互作用在损伤累积方向能诱发一级相变, 但在恢复路径上则维持了连续的二级相变特征, 从而消除了滞回现象.

当同时引入成对与高阶相互作用时, 高阶耦合展现出对滞回行为的显著抑制效应. 当成对耦合强度固定为$ {k}_{\text{p}}=0.75 $, 随着高阶耦合强度$ {k}_{\text{h}} $从0增至0.15, 滞回环消失(如图3(d), (e)所示). 当$ {k}_{\text{p}}= 0.85 $, 增加$ {k}_{\text{h}} $时(图3(g)—(i))也发现了滞回区域宽度有收缩的趋势. 在特定参数组合下(图3(f), $ {k}_{\text{p}}=0.35 $, $ {k}_{\text{h}}=0.3 $), 滞回现象完全消失. 而且, 高阶相互作用的引入还可能导致系统在转变点附近出现复杂的过渡态, 展现出非平凡的中间行为, 对应于稳定的振荡-静息共存的过渡态, 丰富了系统的动力学.

为深入理解其微观机制, 进一步展示了典型状态点对应的时空演化图(图4). 例如, 图4(a)在双爆炸转变曲线上的4个标记点, 对应的参数分别是正向扫描时的$ {P}_{\text{a}}=0.71 $, $ {P}_{\text{a}}=0.83 $以及反向扫描时的$ {P}_{\text{a}}=0.91 $, $ {P}_{\text{a}}=0.83 $. 正向扫描在$ {P}_{\text{a}}=0.71 $时系统为全局同步振荡(图4(c)), 而在正向的$ {P}_{\text{a}}=0.83 $时仍保持全局同步振荡(图4(d)). 反向扫描时$ {P}_{\text{a}}=0.71 $, $ {P}_{\text{a}}=0.83 $对应的则是全局死亡的静息态(图4(e), (f)). 对于无滞回的爆炸性集体失活转变(图4(b)), 其全局振荡与全局死亡的时空演化图分别与图4(c), (e)类似, 为节约正文篇幅, 此处不再重复展示. $ {k}_{\text{p}}=0.5 $, $ {k}_{\text{h}}=0.15 $(图4(g))时转变点附近的过渡态的时空图出现了介于同步振荡与完全死亡之间的部分振荡斑图(图4(h)). 为更精细地刻画这些状态的动力学特征, 图5给出了代表性神经元在3种典型状态下的时间序列, 图中 i 为选取的代表性神经元的编号, 用于展示该神经元的典型时间序列行为. 图5(a)在完全振荡态, 神经元表现出规则、大幅度的周期振荡; 图5(b)在完全振荡死亡态, 膜电位则稳定在静息态; 图5(c)在部分振荡(过渡态), 神经元活动呈现不规则或衰减的振荡模式. 这些时间序列与图4中的时空斑图相互印证, 完整地揭示了系统在不同损伤水平下, 从集体同步振荡、经部分失活的中间态、到最终全局静息的完整动力学演化路径.

为了全面揭示系统老化转变行为对关键控制参数的依赖关系, 并进一步考察高阶耦合强度作为关键调制因子的作用, 我们系统地在损伤比例$ {P}_{\text{a}} $与成对耦合强度$ {k}_{\text{p}} $构成的二维参数平面上进行扫描. 图6—图8分别展示了在$ {k}_{\text{h}}=0 $, $ {k}_{\text{h}}=0.15 $和$ {k}_{\text{h}}=0.3 $时, 系统的序参量$ {P}_{\text{N}} $和$ \Delta {P}_{\text{N}} $在该二维参数空间中的完整响应, 从而得以量化滞回行为并揭示高阶相互作用的效应.

在图6成对耦合中, 序参量滞回差值$ \Delta {P}_{\text{N}} $揭示出两个性质不同的滞回区域. 在较小的$ {k}_{\text{p}} $区域(约$ {k}_{\text{p}} \lt 0.52 $), $ \Delta {P}_{\text{N}} $呈现出一个小于1的非零值区域. 这对应着不彻底的爆炸性转变: 正向扫描中, 系统从全振荡态跃迁至一个部分振荡态, 而非直接进入完全死亡态; 而在反向扫描中, 系统则表现出滞后性, 从而产生滞回现象. 在较大$ {k}_{\text{p}} $区域(约$ {k}_{\text{p}} \gt 0.72 $), $ \Delta {P}_{\text{N}} $出现一个幅值接近1的显著高值带. 这标志着彻底的双爆炸性转变: 正向与反向扫描均发生突然的状态跃迁, 且两条路径分离, 形成滞回环. 同时, 该高值带随$ {k}_{\text{p}} $增大明显向更小的$ {P}_{\text{a}} $方向移动, 表明更强的成对耦合显著降低了系统发生老化转变的损伤阈值, 即网络鲁棒性下降. 值得注意的是, 在中等$ {k}_{\text{p}} $区间(约$ 0.52 \lt {k}_{\text{p}} \lt 0.72 $), $ {P}_{\text{a}} $变化时都有$ \Delta {P}_{N}\approx 0 $, 意味着尽管系统仍可能发生爆炸性失活, 但正向与反向路径基本重合, 滞回现象消失. 这揭示出成对耦合强度本身存在一个临界窗口, 在该窗口内系统可实现无滞回的爆炸性相变.

随着高阶耦合的引入(图7, $ {k}_{\text{h}}=0.15 $), $ \Delta {P}_{\text{N}} $的高值区域在参数空间中明显收缩. 当高阶耦合进一步增强(图8, $ {k}_{\text{h}}=0.3 $), $ \Delta {P}_{\text{N}} $的高值区域进一步缩小. 对比图6—图8可见, 随着高阶耦合强度的增大, 系统对成对耦合强度的依赖性逐渐减弱. 当$ {k}_{\text{h}}=0 $时, $ {P}_{\text{N}} $与$ \Delta {P}_{\text{N}} $的分布均强烈依赖于$ {k}_{\text{p}} $, 滞回区呈明显斜向条带; 当$ {k}_{\text{h}}=0.3 $时, 滞回区几乎消失, $ \Delta {P}_{\text{N}} $的分布趋近于竖直条带, 表明此时系统的宏观状态主要由损伤比例$ {P}_{\text{a}} $决定, 对$ {k}_{\text{p}} $的依赖显著减弱. 这一递变趋势揭示了高阶相互作用在调控网络动力学韧性中的关键作用: 非线性三体协同项提供了额外驱动机制, 使系统部分摆脱成对耦合的主导.

本文通过构建全局耦合的FitzHugh-Nagumo神经元高阶网络模型, 系统揭示了成对与高阶相互作用在调控老化转变中的协同机制与各自独特作用. 研究发现, 成对相互作用是形成爆炸性老化失活与滞回现象的基础: 随着成对耦合强度的变化, 系统表现出从“不彻底”的单爆炸性转变、无滞回的爆炸性转变到“彻底”的双爆炸性强滞回转变的丰富动力学, 且强耦合会导致临界损伤阈值显著左移, 表明网络鲁棒性对耦合强度具有非单调依赖性.

更为重要的是, 高阶相互作用的引入展现了与成对耦合本质不同的调控能力. 在动力学层面, 高阶耦合不仅能诱发单爆炸性转变, 消除恢复路径上的滞后效应, 还能在成对耦合基础上抑制滞回行为, 表现为滞回区域在参数空间中的系统性收缩. 高阶相互作用在老化转变临界区诱导出过渡态, 形成振荡与静息神经元共存的时空斑图. 这些效应共同表明高阶相互作用能增强系统在损伤下的功能韧性: 既通过抑制滞回提升状态恢复的可逆性, 又通过过渡态提供功能缓冲, 避免系统在损伤积累过程中发生急剧转变.

然而, 真实神经系统往往具有更为复杂的拓扑结构(如小世界网络、无标度网络)及更为多样的耦合机制(如化学突触的非线性耦合). 将本文揭示的高阶相互作用动力学效应推广至这些更真实的网络设定中, 探讨其与拓扑异质性的耦合效应、与化学突触非线性的协同机制, 将是未来极具价值的研究方向.