Kuromoto模型

Kuromoto模型(或者Kuramoto–Daido模型), 最早是由蔵本 由紀提出的。 这是一个描述同步现象的模型。更具体的，这是一种描述大量耦合在一起的振子产生同步震荡现象的模型。 最早的模型是为了描述化学与生物学过程中的震荡现象而提出的，后来该模型又被用到了包括神经科学、 火焰震荡动力学等更广泛的领域。 Kuramoto自己也诧异于该模型还能够用来描述一些物理现象，例如一组约瑟夫森结的集体运动。

该模型有很多的假设，包括弱耦合性，振子的同一或弱同一性，每一对振子之间的交互强度依赖于它们在相位差上的正弦值等。

定义

在Kuramoto模型最流行的版本中，每个振荡器都被认为具有其自身的固有自然频率[math]\displaystyle{ \omega_i }[/math]，并且每个振荡器都与其他所有振荡器以相同的方式耦合。惊人的是，这个完全非线性系统的模型在振荡器数量趋近于无穷大（[math]N\rightarrow \infty[/math]）时可以精确求解，或者通过自洽性论证，可以得到极大值参数的稳态解。 该模型最流行的形式是如下形式：

[math]\displaystyle{ \frac{d \theta_i}{d t} = \omega_i + \frac{1}{N} \sum_{j=1}^{N} K_{ij} \sin(\theta_j - \theta_i), \qquad i = 1 \ldots N }[/math]，

其中系统由N个极限环振荡器组成，各自具有相位[math]\displaystyle{ \theta_i }[/math]和耦合常数K。 我们可以在系统中添加噪声. 在这种情况下，原方程修改为：

[math]\displaystyle{ \frac{d \theta_i}{d t} = \omega_{i}+\zeta_{i}+\dfrac{K}{N}\sum_{j=1}^N\sin(\theta_{j}-\theta_{i}) }[/math]

，其中[math]\displaystyle{ \zeta_{i} }[/math]是涨落噪声，且是时间的函数。如果将噪声视为白噪声，则我们有：

[math]\displaystyle{ \langle\zeta_{i}(t)\rangle=0 }[/math] ，

[math]\displaystyle{ \langle\zeta_{i}(t)\zeta_{j}(t')\rangle=2D\delta_{ij}\delta(t-t') }[/math]

其中[math]\displaystyle{ D }[/math]表示噪声的强度。

变换

为了求得该模型（至少在振荡器数量趋于无穷大时）的精确求解，我们可以实施如下变换：

定义有序参数r和ψ，分别表示相位一致性和平均相位，如下所示：

[math]\displaystyle{ re^{i \psi} = \frac{1}{N} \sum_{j=1}^{N} e^{i \theta_j} }[/math]。

这里r表示振荡器群的相位一致性（phase coherence)，ψ表示平均相位。将该变换代入方程后，得到：

[math]\displaystyle{ \frac{d \theta_i}{d t} = \omega_i + K r \sin(\psi-\theta_i) }[/math]。

这样，振荡器的方程不再显式耦合；相反，序参数制约着集体行为。通常，我们会进行进一步的变换：引入一个旋转坐标系，其中所有振荡器的相位统计平均值为零（即[math]\displaystyle{ \psi=0 }[/math]）。最后，方程变为：

[math]\displaystyle{ \frac{d \theta_i}{d t} = \omega_i - K r \sin(\theta_i) }[/math]。

N很大时的极限情况

下面我们考虑当振荡器数量趋于无穷大时的情况。假设固有自然频率的分布为g(ω（已归一化）。然后假设在给定相位θ和给定自然频率ω时，振荡器的密度在时间t为[math]\displaystyle{ \rho(\theta, \omega, t) }[/math]。由于归一化条件，我们有：

[math]\displaystyle{ \int_{-\pi}^{\pi} \rho(\theta, \omega, t) , d \theta = 1. }[/math]

因此，振荡器密度的连续性方程为

[math]\displaystyle{ \frac{\partial \rho}{\partial t} + \frac{\partial}{\partial \theta}[\rho v] = 0, }[/math]

其中v是振荡器的漂移速度，通过在变换后的控制方程中取[math]N\rightarrow \infty[/math]极限得到，使得

[math]\displaystyle{ \frac{\partial \rho}{\partial t} + \frac{\partial}{\partial \theta}[\rho \omega + \rho K r \sin(\psi-\theta)] = 0. }[/math]

最后，必须将有序参数的定义改写为连续（无穷N）极限。[math]\displaystyle{ \theta_i }[/math]必须替换为其在所有[math]\displaystyle{ \omega }[/math]上的系综平均值，并用积分替换求和，得到

[math]\displaystyle{ r e^{i \psi} = \int_{-\pi}^{\pi} e^{i \theta} \int_{-\infty}^{\infty} \rho(\theta, \omega, t) g(\omega) , d \omega , d \theta. }[/math]

大N极限下的解

当所有振荡器以随机漂移的方式呈现不相干（incoherent）状态时，对应的解为[math]\displaystyle{ \rho = 1/(2\pi) }[/math]。此时[math]\displaystyle{ r = 0 }[/math]，振荡器之间没有相位一致性（coherence）。它们均匀地分布在所有可能的相位中，整体处于统计稳态（尽管单个振荡器仍根据其固有ω而改变相位）。

当耦合强度K足够大时，模型可能会出现全同步（fully synchronized）的解。在全同步状态下，所有振荡器共享相同的频率，尽管它们的相位可能不同。

部分同步的解对应的状态是只有部分振荡器（靠近集合平均自然频率的振荡器）实现同步，其余振荡器则以不相干的方式漂移。从数学上来说，这种状态下，对于”锁定“（locked）振荡器的密度为： [math]\displaystyle{ \rho = \delta\left(\theta - \psi - \arcsin\left(\frac{\omega}{K r}\right)\right) }[/math]，

而漂移振荡器的密度为： [math]\displaystyle{ \rho = \frac{\text{归一化常数}}{\omega - K r \sin(\theta - \psi)}。 }[/math] 阶段条件为[math]\displaystyle{ |\omega| \lt K r }[/math]。

当g（自然频率分布）为单峰且对称时，系统的稳定状态解为： [math]\displaystyle{ r=rK \int_{-\pi / 2}^{\pi / 2} \cos ^2 \theta g(K r \sin \theta) d \theta。 }[/math]

随着耦合强度K的增加，当K达到临界值K_c = 2/\pi g(0)时，系统的长期平均[math]\displaystyle{ r=0 }[/math]；但当K稍大于临界值（即K= K_c(1+\mu)，其中μ > 0是一个小量）时，模型的解近似为： [math]\displaystyle{ r \approx \sqrt{\frac{16}{\pi K_{\mathrm{c}}^3}} \sqrt{\frac{\mu}{-g^{\prime \prime}(0)}}。 }[/math]

小N情况

当N较小时，上述解失效，因为无法使用连续近似。

对于N=2的情况是平凡的。在旋转框架中，[math]\displaystyle{ \omega_1 = -\omega_2 }[/math]，因此系统的描述完全由两振荡器之间的角度差[math]\displaystyle{ \Delta \theta = \theta_1 - \theta_2 }[/math]给出。当[math]\displaystyle{ K \lt K_c = 2|\omega_1| }[/math]时，角度会在圆周上循环（即快振荡器会不断绕慢振荡器环行）。当[math]\displaystyle{ K \gt K_c }[/math]时，角度会进入一个稳定的吸引子（即两振荡器的相位锁定）。类似地，N=3的情况下，状态空间是二维环面，系统在二维环面上演化，演化过程无法出现混沌。

混沌最早出现在N=4的情况下。在某些[math]\displaystyle{ \omega_1, \omega_2, \omega_3, K }[/math]的设置下，系统存在奇异吸引子。

与哈密顿系统的联系

耗散型Kuromoto模型包含在某些具有以下形式的哈密顿函数的保守哈密顿系统中：

[math]\displaystyle{ \mathcal{H}(q_1,\ldots, q_N, p_1, \ldots, p_N) = \sum_{i = 1}^N \frac{\omega_i}{2} (q_i^2 + p_i^2) + \frac{K}{4N} \sum_{i,j = 1}^{N} (q_i p_j - q_j p_i) ( q_j^2 + p_j^2 - q_i^2 - p_i^2 ) }[/math]

在将正则变换应用于作用-角度变量（以作用[math]\displaystyle{ I_i = \left( q_i^2 + p_i^2 \right)/2 }[/math]和角度（相位）[math]\displaystyle{ \phi_i = \mathrm{arctan} \left( q_i/p_i \right) }[/math]为变量）后，Kuromoto动力学在常数[math]\displaystyle{ I_i \equiv I }[/math]的不变流形上出现。经过变换的哈密顿函数为：

[math]\displaystyle{ \mathcal{H'}(I_1, \ldots I_N, \phi_1 \ldots, \phi_N) = \sum_{i = 1}^N \omega_i I_i - \frac{K}{N} \sum_{i = 1}^{N} \sum_{j = 1}^{N} \sqrt{I_j I_i} (I_j - I_i) \sin(\phi_j - \phi_i), }[/math]

哈密顿运动方程变为：

[math]\displaystyle{ \frac{d I_i}{dt} = - \frac{\partial \mathcal{H}'}{\partial \phi_i} = - \frac{2 K}{N} \sum_{k=1}^N \sqrt{I_k I_i} (I_k - I_i) \cos(\phi_k - \phi_i) }[/math]

和

[math]\displaystyle{ \frac{d \phi_i}{dt} = \frac{\partial \mathcal{H}'}{\partial I_i} = \omega_i + \frac{K}{N} \sum_{k=1}^N \left[ 2 \sqrt{I_i I_k} \sin(\phi_k - \phi_i) \right. \left. + \sqrt{I_k/I_i} (I_k - I_i) \sin(\phi_k - \phi_i) \right] . }[/math]

因此，当[math]\displaystyle{ I_j = I }[/math]时，流形是守恒的，因为[math]\displaystyle{ \frac{d I_i}{dt} = 0 }[/math]，而相位动力学[math]\displaystyle{ \frac{d \phi_i}{dt} }[/math]成为Kuromoto模型的动力学（当[math]\displaystyle{ I = 1/2 }[/math]时，耦合常数相同）。这一类哈密顿系统包括某些量子-经典系统，例如玻色-爱因斯坦凝聚态。

模型的变体

Kuromoto模型存在多种类型的变体。有些模型改变了拓扑结构，有些允许使用异构权重，而其他一些变体则受到原始模型的启发，但功能形式可能都不同。

在具有不同相位相互作用函数和空间耦合拓扑结构的二维类似Kuramoto振荡器阵列中，存在不同的同步模式：(A) 风车模式，(B) 波模式，(C) 嵌合体模式，(D) 嵌合体与波模式的结合。颜色标度表示振荡器的相位。

网络拓扑结构的变体

除了原来的全连接（all-to-all）拓扑结构外，类似于复杂网络的足够密集的拓扑结构也可以使用原模型解中所用的平均场方法（见上文中的#变换和#大N极限部分了解更多信息）。环形和耦合群体等网络拓扑结构支持奇美拉态。此外，还可以研究在一维拓扑结构（如链和环）中模型的行为。在这种拓扑结构中，耦合无法根据1/N进行缩放，因此无法应用标准的平均场方法，因此必须依赖针对每种情况的分析，尽可能利用对称性，这可能为抽象出一般的解原则提供依据。

在具有扩散局部耦合的二维Kuromoto网络中，我们可以轻松观察到均匀同步、波和螺旋的构型。使用Turings稳定性分析的方法可以分析确定这些波的稳定性。当局部耦合在所有地方为正时，均匀同步往往是稳定的；而当负的长程连接（抑制性环绕耦合）存在时，则会产生波。波和同步通过一种名为纹波（ripple）的不同拓扑分支的解决方案相连。这些是从均匀状态（或波状态）通过霍普夫二倍分支出现的低幅度空间周期偏差。Wiley、Strogatz和Michelle Girvan预测了纹波解的存在（但未观察到），并称其为多重曲线q态。

可以通过使用适应性模型使Kuromoto模型的拓扑结构能够自适应调整，从而以自组织方式增强同步和连通性。

具有最小度数至少为[math]\displaystyle{ d_{min}\ge 0.5\ n }[/math]的图将是连通的，然而，为了使图同步得更好，需要满足一个临界连接度阈值[math]\displaystyle{ \mu_c }[/math]，使得任何具有[math]\displaystyle{ n }[/math]个节点且最小度数[math]\displaystyle{ d_{min}\ge \mu_c (n-1) }[/math]的图必须实现全局同步（当[math]\displaystyle{ n }[/math]足够大时）。已知[math]\displaystyle{ \mu_c }[/math]的最小值和最大值介于[math]\displaystyle{ 0.6875\le \mu_c\le 0.75 }[/math]之间。

类似地，已知ER随机图在节点数趋近于无穷大时，边概率为[math]\displaystyle{ p=(1+\epsilon)\ln (n)/n }[/math]时将是连通的，并且已经提出猜想认为此值也是这些随机图实现同步的节点数，2022年的一篇预印本声称已证明了这一点。

网络拓扑和网络权重的变体：从车辆协调到脑部同步

控制领域的一些研究集中于具有异构权重的网络中的Kuromoto模型（即任意两个振荡器之间的互连强度可以是任意值）。此模型的动力学如下所示:

[math]\displaystyle{ \frac{d \theta_i}{d t} = \omega_i + \sum_{j=1}^{N} a_{ij} \sin(\theta_j - \theta_i), \qquad i = 1 \ldots N }[/math]，

其中[math]\displaystyle{ a_{ij} }[/math]是非零正实数，当振荡器[math]\displaystyle{ j }[/math]连接到振荡器[math]\displaystyle{ i }[/math]时。这种模型允许对例如电网、成群行为、鱼群和车辆协调等更为现实的现象进行研究。在Dörfler及其同事的工作中，他们提出了多个定理提供了严格的条件，以实现此模型的相位和频率同步。此外，受神经科学实验观察的启发，进一步的研究侧重于在任意网络拓扑上推导异构Kuromoto振荡器的集群同步的分析条件。由于Kuromoto模型似乎在评估大脑中的同步现象方面起着关键作用，支持实证结果的理论条件可能有助于深入理解神经同步现象。

相位交互函数的变体

Kuromoto第一个使用傅里叶分量方法对对两个振荡器之间的相位相互作用进行了近似求解，即[math]\displaystyle{ \Gamma(\phi) = \sin(\phi) }[/math]，其中[math]\displaystyle{ \phi = \theta_j - \theta_i }[/math]。通过包含更高阶的傅里叶分量，我们可以获得更好的近似：

[math]\displaystyle{ \Gamma(\phi) = \sin(\phi) + a_1 \sin(2\phi + b_1) + ... + a_n \sin(2n\phi + b_n) }[/math]，

其中参数[math]\displaystyle{ a_i }[/math]和[math]\displaystyle{ b_i }[/math]必须通过估计得到。例如，可以使用保留交互函数的前四个傅里叶分量的耦合振荡器来复制一组弱耦合的Hodgkin–Huxley神经元的同步。引入更高阶的相位交互项也会引发有趣的动态现象，例如部分同步状态、异离心循环，以及混沌动力学。

可用性

pyclustering 库包含了库拉莫托模型及其修改版本的 Python 和 C++ 实现。此外，该库还包含基于库拉莫托模型和相位振荡器的振荡网络（用于集群分析、模式识别、图着色和图像分割）。