Hopf模型

Hopf模型（Hopf model）是以 Hopf 分岔正规形为基础的低维振子模型，用来描述系统从衰减振荡过渡到持续振荡的过程。正规形指在分岔点附近，通过变量变换提取关键动力学后得到的简化方程。Hopf分岔是指平衡点的一对共轭复特征值穿过虚轴时，附近产生或消失周期轨道的局部分岔。本词条以常用的超临界Hopf模型为主：控制参数超过阈值后，系统出现振幅从零连续增大的稳定极限环，即能够吸引邻近轨迹的周期运动。

在神经动力学中，Hopf 模型既可描述神经系统接近振荡起始点时的局部行为，也可作为全脑网络的区域动力学单元。每个区域保留振幅和相位两个自由度，再通过结构连接与其他区域耦合。用于静息态功能磁共振成像时，一种常用版本直接以模型的实部表示血氧水平依赖信号（BOLD）的低频波动，属于对观测信号动力学的简化描述。

基本信息

项目	内容
英文名称	Hopf model；常用三次正规形也称 Stuart–Landau 振子
模型类型	连续时间非线性振子，可加入外部输入、噪声和网络耦合
状态变量	一个复变量 [math]\displaystyle{ z=x+\mathrm{i}y }[/math]，等价于两个实变量 [math]\displaystyle{ x,y }[/math]
主要参数	分岔参数 [math]\displaystyle{ a }[/math]、本征角频率 [math]\displaystyle{ \omega }[/math]、振幅饱和系数 [math]\displaystyle{ b }[/math]
典型行为	稳定焦点、衰减振荡、稳定极限环、共振、频率锁定
常见用途	振荡起始机制、神经系统共振、结构连接约束下的全脑活动建模

理论渊源与发展

Hopf 模型来自非线性微分方程的局部分岔理论。当系统在平衡点附近只有一对关键特征值接近虚轴，而其余方向保持稳定时，可以把长期局部行为约化到二维中心流形，即容纳这些关键运动方向的不变曲面。正规形保留振荡的产生、振幅增长和非线性饱和机制，使不同高维系统可以用相似的低维方程分析。其有效范围主要是分岔点附近。

这一思想在计算神经科学中用于理解神经元的振荡与共振。Izhikevich 在 2000 年的工作中分析了 Hopf 分岔附近神经系统对输入频率和刺激间隔的选择性响应，并将这类行为与神经元的兴奋性联系起来。在全脑建模中，Deco 等人于 2017 年将带噪声的 Hopf 振子放在经验结构连接网络上，同时比较静态功能连接、功能连接的时间变化和相位同步的波动。

单个振子的数学形式

复数形式与实数形式

一个无外部输入的超临界 Hopf 振子可写为：

[math]\displaystyle{ \frac{\mathrm{d}z}{\mathrm{d}t}=(a+\mathrm{i}\omega)z-b|z|^2z,\qquad b\gt 0,\quad\omega\ne0. }[/math]

(1)

其中 [math]\displaystyle{ \mathrm{i}^2=-1 }[/math]，[math]\displaystyle{ z=x+\mathrm{i}y }[/math]，[math]\displaystyle{ |z|=\sqrt{x^2+y^2} }[/math]。参数 [math]\displaystyle{ a }[/math] 决定小振幅扰动增长还是衰减，[math]\displaystyle{ \omega }[/math] 决定旋转速度，[math]\displaystyle{ b }[/math] 使振幅在增大后受到抑制。若时间以秒计，角频率的单位为弧度每秒，普通频率为 [math]\displaystyle{ f=|\omega|/(2\pi) }[/math] Hz。全脑模型常通过状态尺度选择把 [math]\displaystyle{ b }[/math] 设为 1。

把式（1）的实部和虚部分开，得到：

[math]\displaystyle{ \begin{aligned}\dot{x}&=[a-b(x^2+y^2)]x-\omega y,\\\dot{y}&=[a-b(x^2+y^2)]y+\omega x.\end{aligned} }[/math]

(2)

此处 [math]\displaystyle{ x,y }[/math] 是振子的两个状态分量。它们的物理含义由建模层级和观测方式决定，例如全脑 BOLD 模型把 [math]\displaystyle{ x }[/math] 作为与区域 BOLD 波动比较的信号。用于更复杂神经系统的局部约化时，状态变量还需通过原系统的变量变换解释。

振幅、相位与稳定性

令 [math]\displaystyle{ z=r\exp(\mathrm{i}\theta) }[/math]，其中 [math]\displaystyle{ r\ge0 }[/math] 是振幅，[math]\displaystyle{ \theta }[/math] 是相位。对 [math]\displaystyle{ r\gt 0 }[/math]，式（1）给出：

[math]\displaystyle{ \dot{r}=ar-br^3,\qquad\dot{\theta}=\omega. }[/math]

(3)

式（3）把问题分为径向增长和相位旋转。雅可比矩阵由方程的一阶偏导组成，用于描述小扰动的线性演化。原点处的矩阵及其特征值为：

[math]\displaystyle{ \mathbf{J}_0=\begin{pmatrix}a&-\omega\\\omega&a\end{pmatrix},\qquad\lambda_{\pm}=a\pm\mathrm{i}\omega. }[/math]

(4)

• [math]\displaystyle{ a\lt 0 }[/math]：原点是稳定焦点，轨迹一边旋转一边靠近原点。小振幅的衰减近似为 [math]\displaystyle{ r(t)\approx r(0)e^{at} }[/math]。
• [math]\displaystyle{ a\gt 0 }[/math]：原点失稳，所有非零初始振幅趋向稳定的非零振幅，形成周期运动。精确位于原点的无噪声轨迹仍停留在原点。
• [math]\displaystyle{ a=0 }[/math]：线性特征值实部为零，稳定性要由非线性项判断。在式（1）中，三次项仍使振幅衰减，原点保持非线性渐近稳定。

对 [math]\displaystyle{ a\gt 0 }[/math]，极限环振幅和周期为：

[math]\displaystyle{ r_*=\sqrt{\frac{a}{b}},\qquad T=\frac{2\pi}{|\omega|}. }[/math]

(5)

振幅在阈值附近按 [math]\displaystyle{ \sqrt{a} }[/math] 增长，频率在此简化模型中保持不变。若在式（1）的三次项中把 [math]\displaystyle{ b }[/math] 扩展为 [math]\displaystyle{ b+\mathrm{i}c }[/math]，振幅方程保持原形式，相位速度变为 [math]\displaystyle{ \dot\theta=\omega-cr^2 }[/math]，从而允许频率随振幅变化。正负 [math]\displaystyle{ a }[/math] 的简单判据属于所选单振子正规形，分析复杂网络时还需检查整个网络的线性化。

在分岔点 [math]\displaystyle{ a=0 }[/math]，对初始振幅 [math]\displaystyle{ r_0\gt 0 }[/math]，式（3）可直接积分为：

[math]\displaystyle{ r(t)=\frac{r_0}{\sqrt{1+2br_0^2t}}. }[/math]

(6)

这里的衰减呈幂律，比稳定焦点中的指数衰减缓慢。靠近阈值时，[math]\displaystyle{ a\lt 0 }[/math] 的小振幅弛豫时间约为 [math]\displaystyle{ 1/|a| }[/math]，即扰动衰减所需的特征时间。[math]\displaystyle{ a\gt 0 }[/math] 时，极限环附近的径向弛豫时间约为 [math]\displaystyle{ 1/(2a) }[/math]。这说明接近分岔点会使恢复过程变慢。

超临界与次临界

式（1）的负三次项产生超临界分岔。相反，若三次项为正，在阈值的稳定平衡点一侧可以出现不稳定的小周期轨道，对应次临界 Hopf 分岔。要描述次临界模型中较大振幅的稳定振荡，通常还需加入使振幅饱和的更高阶项。跳变与滞后是否出现取决于这些额外项及整体分岔结构。

一般系统的 Hopf 分岔判定还要求共轭特征值以非零速度穿过虚轴，以及控制三次非线性的系数非零。这些条件保证分岔的通常形式成立。单纯观察到一段振荡时间序列，还不足以确定其分岔类型。

数值示例一：从衰减振荡到极限环

图1使用式（2），取 [math]\displaystyle{ b=1,\omega=1 }[/math]，初始条件 [math]\displaystyle{ (x_0,y_0)=(0.8,0) }[/math]，比较 [math]\displaystyle{ a=-0.2,0,0.2 }[/math]。时间为无量纲时间，计算区间为 [math]\displaystyle{ 0\le t\le40 }[/math]。

图1：超临界 Hopf 模型的分岔与轨迹。（A）长期振幅随控制参数变化。灰实线是稳定原点，灰虚线是不稳定原点，橙线是稳定周期轨道的振幅。（B）三组参数的振幅随时间变化：蓝实线 [math]\displaystyle{ a=-0.2 }[/math]，绿点线 [math]\displaystyle{ a=0 }[/math]，橙虚线 [math]\displaystyle{ a=0.2 }[/math]。（C）负参数下轨迹向原点螺旋收敛。（D）正参数下轨迹趋向半径约 0.4472 的极限环，灰虚线标出该圆。相图的空心点表示起点，实心点表示 [math]\displaystyle{ t=40 }[/math] 的终点。[math]\displaystyle{ x,y,r }[/math] 及时间均为无量纲量。来源：杨明哲依据式（2）进行数值计算并绘制，CC BY-SA 3.0。

数值结果中，[math]\displaystyle{ a=-0.2 }[/math] 的振幅趋近零，[math]\displaystyle{ a=0.2 }[/math] 的振幅趋近 [math]\displaystyle{ \sqrt{0.2}\approx0.4472 }[/math]。[math]\displaystyle{ a=0 }[/math] 时仍可在有限时间内观察到旋转，但振幅按式（6）逐渐衰减。因而分岔点处的暂态旋转与稳定的自持振荡应结合长期振幅来区分。

外部输入、共振与噪声

在式（1）右侧加入 [math]\displaystyle{ u(t)=F\exp(\mathrm{i}\Omega t) }[/math]，其中 [math]\displaystyle{ F\ge0 }[/math] 是输入强度，[math]\displaystyle{ \Omega }[/math] 是输入角频率。对稳定焦点 [math]\displaystyle{ a\lt 0 }[/math]，当稳态响应足够小、三次项可以忽略时，响应振幅近似为：

[math]\displaystyle{ A\approx\frac{F}{\sqrt{a^2+(\Omega-\omega)^2}}. }[/math]

(7)

输入频率接近本征频率时，响应更强，这就是共振。式（7）适用于 [math]\displaystyle{ bA^2\ll|a| }[/math] 的线性响应区。在 [math]\displaystyle{ a=0 }[/math] 且 [math]\displaystyle{ \Omega=\omega }[/math] 时，三次项决定稳态振幅，得到 [math]\displaystyle{ A=(F/b)^{1/3} }[/math]。因此临界点的响应需要保留非线性饱和项。神经系统接近 Hopf 分岔时的频率选择性可以用这类机制理解。

加入随机驱动后，即使确定性系统在 [math]\displaystyle{ a\lt 0 }[/math] 时最终回到原点，持续噪声也能反复激起衰减振荡。这样的噪声驱动振荡由外部随机输入维持，其振幅和相位不断波动。自持极限环则能够在无噪声条件下维持有限振幅。两种情况都可以出现谱峰，需要结合参数、动力学和数据共同解释。

耦合网络与全脑模型

网络方程

令 [math]\displaystyle{ z_i }[/math] 表示第 [math]\displaystyle{ i }[/math] 个节点，带扩散耦合和独立加性噪声的网络可写为：

[math]\displaystyle{ \mathrm{d}z_i=\left[(a_i+\mathrm{i}\omega_i)z_i-b_i|z_i|^2z_i+G\sum_{j=1}^{N}C_{ij}(z_j-z_i)\right]\mathrm{d}t+\sigma_i(\mathrm{d}W_{x,i}+\mathrm{i}\,\mathrm{d}W_{y,i}). }[/math]

(8)

其中 [math]\displaystyle{ N }[/math] 是节点数，[math]\displaystyle{ \mathbf{C}=(C_{ij}) }[/math] 是非负连接权重矩阵，约定 [math]\displaystyle{ C_{ij} }[/math] 表示节点 [math]\displaystyle{ j }[/math] 对节点 [math]\displaystyle{ i }[/math] 的连接，[math]\displaystyle{ C_{ii}=0 }[/math]。[math]\displaystyle{ G\ge0 }[/math] 是全局耦合强度，[math]\displaystyle{ b_i\gt 0 }[/math]。[math]\displaystyle{ W_{x,i},W_{y,i} }[/math] 是彼此独立的标准维纳过程，即增量方差等于时间间隔的随机过程，[math]\displaystyle{ \sigma_i }[/math] 控制噪声强度。式（8）明确采用白噪声的随机微分方程约定，数值实现见式（13）。

扩散项根据节点状态的差异产生相互作用。将式（8）的确定性部分换为极坐标，可以得到相位方程：

[math]\displaystyle{ \dot\theta_i=\omega_i+G\sum_j C_{ij}\frac{r_j}{r_i}\sin(\theta_j-\theta_i),\qquad r_i\gt 0. }[/math]

(9)

当振幅近似相等并保持稳定时，相位耦合简化为相位差的正弦函数。此时可以联系仅保留相位的 Kuramoto 模型，即描述不同频率振子如何协调相位的模型。接近振荡阈值时，振幅变化较大，Hopf 模型保留的振幅自由度尤其有用。存在噪声时，相位还会受到随机驱动。

局部分岔与网络稳定性

节点间耦合会改变集体模态，网络的稳定性需要由整体雅可比矩阵判断。对参数同质、连接矩阵对称且权重非负的无噪声网络，定义拉普拉斯矩阵 [math]\displaystyle{ \mathbf{L}=\mathbf{D}-\mathbf{C} }[/math]，其中 [math]\displaystyle{ D_{ii}=\sum_jC_{ij} }[/math]，其余元素为零。若 [math]\displaystyle{ \ell_k }[/math] 是 [math]\displaystyle{ \mathbf{L} }[/math] 的特征值，原点附近各模态的特征值为：

[math]\displaystyle{ \lambda_{k,\pm}=a-G\ell_k\pm\mathrm{i}\omega. }[/math]

(10)

因为 [math]\displaystyle{ \ell_k\ge0 }[/math]，扩散耦合在这里抑制空间差异模态，而 [math]\displaystyle{ \ell_0=0 }[/math] 的同步模态仍由 [math]\displaystyle{ a }[/math] 控制。因此，在这个特定同质网络中，增大 [math]\displaystyle{ G }[/math] 本身不会使所有 [math]\displaystyle{ a\lt 0 }[/math] 的节点产生确定性自持振荡。异质频率、不同局部参数和其他耦合形式需要分别分析。

数值示例二：两个振子的频率锁定

考虑式（8）的无噪声双节点网络，取 [math]\displaystyle{ a_1=a_2=b_1=b_2=1 }[/math]，[math]\displaystyle{ \omega_1=1,\omega_2=1.4 }[/math]，[math]\displaystyle{ C_{12}=C_{21}=1 }[/math]。初始振幅均为 1，初始相位为 [math]\displaystyle{ \theta_1=0,\theta_2=0.7 }[/math]。该对称设置下，两振幅在演化中保持相等，相位差 [math]\displaystyle{ \delta=\theta_2-\theta_1 }[/math] 满足：

[math]\displaystyle{ \dot\delta=\Delta\omega-2G\sin\delta,\qquad\Delta\omega=\omega_2-\omega_1. }[/math]

(11)

当 [math]\displaystyle{ |\Delta\omega|\lt 2G }[/math] 时，存在稳定的相位锁定解 [math]\displaystyle{ \delta_*=\arcsin[\Delta\omega/(2G)] }[/math]，其余径向条件也需保证振幅非零。本例取 [math]\displaystyle{ G=0.35 }[/math]，得到 [math]\displaystyle{ \delta_*\approx0.6082 }[/math] 弧度。两个振子的共同角频率为 [math]\displaystyle{ (\omega_1+\omega_2)/2=1.2 }[/math]，共同振幅为 [math]\displaystyle{ \sqrt{1+G(\cos\delta_*-1)}\approx0.9681 }[/math]。

图2：扩散耦合产生频率锁定。（A）[math]\displaystyle{ G=0 }[/math] 时，两振子分别保持角频率 1 和 1.4。（B）[math]\displaystyle{ G=0.35 }[/math] 时，两者以共同角频率 1.2 振荡，并保持非零相位差。蓝实线和橙虚线分别表示振子1和振子2的实部。（C）连续展开后的相位差，即跨越 [math]\displaystyle{ 2\pi }[/math] 后仍累计的相位差。蓝线为无耦合情形，橙虚线为耦合情形，灰点线标出锁定值约 0.6082 弧度。（A）、（B）显示 [math]\displaystyle{ 20\le t\le40 }[/math]，（C）显示完整计算区间。时间及活动幅值为无量纲量。来源：杨明哲依据式（8）进行数值计算并绘制，CC BY-SA 3.0。

无耦合时，相位差按 [math]\displaystyle{ 0.7+0.4t }[/math] 增长。加入耦合后，相位差趋向常数，频率较快的振子减速，较慢的振子加速。频率锁定表示两者长期频率相等，同时可以保留固定的非零相位差。本例的两条波形因此保持错位。该示例用于说明耦合机制，参数采用无量纲设置。

全脑建模的研究案例

静息态活动与亚稳态波动

Deco 等人的 2017 年研究使用 24 名健康参与者的数据，将大脑皮层划分为 68 个区域，以弥散成像得到的结构连接耦合区域 Hopf 振子。模型根据各区域 BOLD 信号的谱峰设置本征频率，比较约 0.04–0.07 Hz 频段的活动。这里的频率 [math]\displaystyle{ f_i }[/math] 与式（8）的角频率满足 [math]\displaystyle{ \omega_i=2\pi f_i }[/math]。例如 0.05 Hz 对应约 20 秒的周期。

该研究同时考察功能连接（FC，区域信号之间的统计相关）与功能连接动态（FCD，不同时间窗的连接模式有多相似）。其 FCD 使用 60 秒时间窗，相邻窗口重叠 40 秒，再比较窗口间 FC 矩阵上三角元素的相关性。静态 FC 衡量整段记录的平均关系，FCD 则能反映连接模式随时间的变化。

研究还用 Kuramoto 序参量衡量相位集中程度，并将其时间标准差称为亚稳态性（metastability），即全局相位同步程度随时间波动的幅度：

[math]\displaystyle{ R(t)=\left|\frac{1}{N}\sum_{j=1}^{N}e^{\mathrm{i}\phi_j(t)}\right|,\qquad M=\operatorname{std}_t[R(t)]. }[/math]

(12)

其中 [math]\displaystyle{ \phi_j(t) }[/math] 是信号相位，[math]\displaystyle{ 0\le R(t)\le1 }[/math]，[math]\displaystyle{ R=1 }[/math] 表示所有相位相同。[math]\displaystyle{ M }[/math] 描述同步的时间变化，在这一研究中用于表征亚稳态性。常用的信号相位由窄带滤波后进行希尔伯特变换提取，即把实信号与其相移分量组合以估计瞬时相位。

在所研究的数据与参数设置下，局部参数接近振荡阈值的网络较好地匹配了经验活动的多项统计特征。其解释是：局部动力学、结构连接与噪声共同影响短暂同步模式的形成和消退。这一结果说明 Hopf 模型可用于检验静息态动力学的候选机制，最优参数仍依赖连接权重尺度、噪声及观测处理。

清醒与慢波睡眠的比较

Jobst 等人同年利用 18 名参与者的清醒与慢波睡眠数据，建立了包含 90 个皮层及皮层下区域的 Hopf 网络。慢波睡眠指非快速眼动睡眠中较深的 N3 阶段。研究同时搜索全局耦合强度与统一的局部分岔参数，比较模型和经验 FC 的差异。

拟合结果显示，清醒态较接近衰减振荡与自持振荡的过渡区域，慢波睡眠态则移向更负的分岔参数，表现为更稳定的焦点动力学。研究还通过优化连接权重发现，睡眠态的模型有效相互作用整体降低。作者据此提出，局部稳定性增加与区域间相互作用减弱共同限制了扰动传播。这里的有效连接是特定模型估计的区域相互作用，需与原始结构连接以及直接的干预测量分别理解。

数值实现与参数解释

本文两幅图中的轨迹使用 DOP853 自适应积分，输出间隔为 0.01，最大积分步长为 0.05，相对容差为 [math]\displaystyle{ 10^{-11} }[/math]，绝对容差为 [math]\displaystyle{ 10^{-12} }[/math]。单振子结果已与振幅解析解比较，双振子结果已与式（11）的锁定解比较。双振子计算将最大步长减半后，所报告结果保持一致。图1（A）的分支由式（5）及原点稳定性直接给出。

对式（8）的白噪声约定，可以用 Euler–Maruyama 方法实现，即在确定性 Euler 更新之外加入方差与步长成正比的随机增量：

[math]\displaystyle{ z_i^{n+1}=z_i^n+hF_i(\boldsymbol{z}^n)+\sigma_i\sqrt{h}(\epsilon_{x,i}^n+\mathrm{i}\epsilon_{y,i}^n). }[/math]

(13)

其中 [math]\displaystyle{ h }[/math] 为时间步长，[math]\displaystyle{ \boldsymbol{z}=(z_1,\ldots,z_N)^{\mathsf T} }[/math]，[math]\displaystyle{ F_i }[/math] 是式（8）方括号中的确定性漂移项，[math]\displaystyle{ \epsilon_{x,i}^n,\epsilon_{y,i}^n }[/math] 是各节点、分量及时间步独立的标准正态随机数。噪声增量按 [math]\displaystyle{ \sqrt h }[/math] 缩放，与每个时间步直接加入固定方差噪声的约定不同。比较不同软件或论文参数时，需要同时检查时间单位和噪声定义。

全脑建模还应记录区域划分、结构连接的归一化、频率估计方法、观测信号定义、滤波和采样设置。改变 [math]\displaystyle{ \mathbf C }[/math] 的整体尺度会改变同一 [math]\displaystyle{ G }[/math] 所代表的有效耦合。例如把连接权重整体乘以常数，同时将 [math]\displaystyle{ G }[/math] 除以该常数，式（8）的耦合项保持不变。对随机模型应报告随机种子、重复模拟及结果变动范围，参数拟合与模型评估宜采用相同的信号处理流程。

与其他神经动力学模型的关系及局限

模型	主要状态与解释	与 Hopf 模型的联系
FitzHugh-Nagumo模型	快变量与恢复变量，描述兴奋性及振荡	在出现 Hopf 分岔的参数区间，可用正规形分析局部振荡起始
Wilson-Cowan模型	兴奋性与抑制性群体的平均活动	群体相互作用可产生 Hopf 分岔，正规形概括其局部周期活动
动力学平均场模型	区域群体活动与突触门控变量，包含更具体的突触和输入机制	都可用于结构连接约束的全脑模拟，参数的生理解释和观测映射随模型而异
Hopf模型	振幅与相位两个自由度	便于分析振荡阈值、频率协调与噪声作用

Hopf 正规形以分岔附近的关键动力学为重点。离分岔点较远时，原系统的高阶项、多种振荡模态和其他状态变量可能变得重要。式（1）的单个无噪声自治振子主要产生平衡态或周期运动，更复杂的尖峰波形、簇发或不规则网络活动需要额外机制及相应分析。

用于 BOLD 建模时，状态变量和参数首先服务于观测信号的动力学描述。将模型参数解释为膜电导、特定突触受体或神经元放电率，需要进一步建立生理映射。对经验 FC 或 FCD 的良好拟合支持该模型作为候选解释，也应结合频谱、时间变化和独立数据评估其适用性。

署名

贡献者：杨明哲。