Kuromoto模型

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

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

定义

文件:KuramotoModelPhaseLocking.ogv

Phase locking in the Kuramoto model

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

[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表示振荡器群的相位共质，ψ表示平均相位。将该变换代入方程后，得到：

[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]。

Large N limit

Now consider the case as N tends to infinity. Take the distribution of intrinsic natural frequencies as g(ω) (assumed normalized). Then assume that the density of oscillators at a given phase θ, with given natural frequency ω, at time t is [math]\displaystyle{ \rho(\theta, \omega, t) }[/math]. Normalization requires that

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

The continuity equation for oscillator density will be

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

where v is the drift velocity of the oscillators given by taking the infinite-N limit in the transformed governing equation, such that

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

Finally, the definition of the order parameters must be rewritten for the continuum (infinite N) limit. [math]\displaystyle{ \theta_i }[/math] must be replaced by its ensemble average (over all [math]\displaystyle{ \omega }[/math]) and the sum must be replaced by an integral, to give

[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]

Solutions for the large N limit

The incoherent state with all oscillators drifting randomly corresponds to the solution [math]\displaystyle{ \rho = 1/(2\pi) }[/math]. In that case [math]\displaystyle{ r = 0 }[/math], and there is no coherence among the oscillators. They are uniformly distributed across all possible phases, and the population is in a statistical steady-state (although individual oscillators continue to change phase in accordance with their intrinsic ω).

When coupling K is sufficiently strong, a fully synchronized solution is possible. In the fully synchronized state, all the oscillators share a common frequency, although their phases can be different.

A solution for the case of partial synchronization yields a state in which only some oscillators (those near the ensemble's mean natural frequency) synchronize; other oscillators drift incoherently. Mathematically, the state has

[math]\displaystyle{ \rho = \delta\left(\theta - \psi - \arcsin\left(\frac{\omega}{K r}\right)\right) }[/math]

for locked oscillators, and

[math]\displaystyle{ \rho = \frac{\rm{normalization \; constant}}{(\omega - K r \sin(\theta - \psi))} }[/math]

for drifting oscillators. The cutoff occurs when [math]\displaystyle{ |\omega| \lt K r }[/math].

When [math]\displaystyle{ g }[/math] is unimodal and symmetric, then a stable state solution for the system is [math]\displaystyle{ r=rK \int_{-\pi / 2}^{\pi / 2} \cos ^2 \theta g(K r \sin \theta) d \theta }[/math]As coupling increases, there is a critical value [math]\displaystyle{ K_c = 2/\pi g(0) }[/math] such that when [math]\displaystyle{ K \lt K_c }[/math], the long-term average of [math]\displaystyle{ r=0 }[/math], but when [math]\displaystyle{ K= K_c(1+\mu) }[/math], where [math]\displaystyle{ \mu \gt 0 }[/math] is small, then [math]\displaystyle{ r \approx \sqrt{\frac{16}{\pi K_{\mathrm{c}}^3}} \sqrt{\frac{\mu}{-g^{\prime \prime}(0)}} }[/math].

Small N cases

When N is small, the solutions given above breaks down, as the continuum approximation cannot be used.

The N=2 case is trivial. In the rotating frame [math]\displaystyle{ \omega_1 = -\omega_2 }[/math], and so the system is described exactly by the angle between the two oscillators: [math]\displaystyle{ \Delta \theta = \theta_1 - \theta_2 }[/math]. When [math]\displaystyle{ K \lt K_c = 2|\omega_1| }[/math], the angle cycles around the circle (that is, the fast oscillator keeps lapping around the slow oscillator). When [math]\displaystyle{ K \gt K_c }[/math], the angle falls into a stable attractor (that is, the two oscillators lock in phase). Similarly, the state space of the N=3 case is a 2-dimensional torus, and so the system evolves as a flow on the 2-torus, which cannot be chaotic.

Chaos first occurs when N=4. For some settings of [math]\displaystyle{ \omega_1, \omega_2, \omega_3, K }[/math], the system has a strange attractor.

Connection to Hamiltonian systems

The dissipative Kuramoto model is contained in certain conservative Hamiltonian systems with Hamiltonian of the form

[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]

After a canonical transformation to action-angle variables with actions [math]\displaystyle{ I_i = \left( q_i^2 + p_i^2 \right)/2 }[/math] and angles (phases) [math]\displaystyle{ \phi_i = \mathrm{arctan} \left( q_i/p_i \right) }[/math], exact Kuramoto dynamics emerges on invariant manifolds of constant [math]\displaystyle{ I_i \equiv I }[/math]. With the transformed Hamiltonian

[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]

Hamilton's equation of motion become

[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]

and

[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]

So the manifold with [math]\displaystyle{ I_j = I }[/math] is invariant because [math]\displaystyle{ \frac{d I_i}{dt} = 0 }[/math] and the phase dynamics [math]\displaystyle{ \frac{d \phi_i}{dt} }[/math] becomes the dynamics of the Kuramoto model (with the same coupling constants for [math]\displaystyle{ I = 1/2 }[/math]). The class of Hamiltonian systems characterizes certain quantum-classical systems including Bose–Einstein condensates.

Variations of the models

文件:Spatial synchronization patterns in variants of the Kuramoto Model.png

Distinct synchronization patterns in a two-dimensional array of Kuramoto-like oscillators with differing phase interaction functions and spatial coupling topologies. (A) Pinwheels. (B) Waves. (C) Chimeras. (D) Chimeras and waves combined. Color scale indicates oscillator phase.

There are a number of types of variations that can be applied to the original model presented above. Some models change the topological structure, others allow for heterogeneous weights, and other changes are more related to models that are inspired by the Kuramoto model but do not have the same functional form.

Variations of network topology

Beside the original model, which has an all-to-all topology, a sufficiently dense complex network-like topology is amenable to the mean-field treatment used in the solution of the original model (see Transformation and Large N limit above for more info). Network topologies such as rings and coupled populations support chimera states. One also may ask for the behavior of models in which there are intrinsically local, like one-dimensional topologies which the chain and the ring are prototypical examples. In such topologies, in which the coupling is not scalable according to 1/N, it is not possible to apply the canonical mean-field approach, so one must rely upon case-by-case analysis, making use of symmetries whenever it is possible, which may give basis for abstraction of general principles of solutions.

Uniform synchrony, waves and spirals can readily be observed in two-dimensional Kuramoto networks with diffusive local coupling. The stability of waves in these models can be determined analytically using the methods of Turing stability analysis. Uniform synchrony tends to be stable when the local coupling is everywhere positive whereas waves arise when the long-range connections are negative (inhibitory surround coupling). Waves and synchrony are connected by a topologically distinct branch of solutions known as ripple. These are low-amplitude spatially-periodic deviations that emerge from the uniform state (or the wave state) via a Hopf bifurcation. The existence of ripple solutions was predicted (but not observed) by Wiley, Strogatz and Girvan, who called them multi-twisted q-states.

The topology on which the Kuramoto model is studied can be made adaptive by use of fitness model showing enhancement of synchronization and percolation in a self-organised way.

A graph with the minimal degree at least [math]\displaystyle{ d_{min}\ge 0.5\ n }[/math] will be connected nevertheless for a graph to synchronize a little more it is required for such case it is known that there is critical connectivity threshold [math]\displaystyle{ \mu_c }[/math] such that any graph on [math]\displaystyle{ n }[/math] nodes with minimum degree [math]\displaystyle{ d_{min}\ge \mu_c (n-1) }[/math] must globally synchronise.for [math]\displaystyle{ n }[/math] large enough. The minimum maximum are known to lie between [math]\displaystyle{ 0.6875\le \mu_c\le 0.75 }[/math].

Similarly it is known that Erdős-Rényi graphs with edge probability precisely [math]\displaystyle{ p=(1+\epsilon)\ln (n)/n }[/math] as [math]\displaystyle{ n }[/math] goes to infinity will be connected and it has been conjectured that this value is too the number at which these random graphs undergo synchronization which a 2022 preprint claims to have proved.

Variations of network topology and network weights: from vehicle coordination to brain synchronization

文件:Wesphysdemo - Synchronized Metronomes.webm

Metronomes, initially out of phase, synchronize through small motions of the base on which they are placed. This system has been shown to be equivalent to the Kuramoto model.

Some works in the control community have focused on the Kuramoto model on networks and with heterogeneous weights (i.e. the interconnection strength between any two oscillators can be arbitrary). The dynamics of this model reads as follows:

[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]

where [math]\displaystyle{ a_{ij} }[/math] is a nonzero positive real number if oscillator [math]\displaystyle{ j }[/math] is connected to oscillator [math]\displaystyle{ i }[/math]. Such model allows for a more realistic study of, e.g., power grids, flocking, schooling, and vehicle coordination. In the work from Dörfler and colleagues, several theorems provide rigorous conditions for phase and frequency synchronization of this model. Further studies, motivated by experimental observations in neuroscience, focus on deriving analytical conditions for cluster synchronization of heterogeneous Kuramoto oscillators on arbitrary network topologies. Since the Kuramoto model seems to play a key role in assessing synchronization phenomena in the brain, theoretical conditions that support empirical findings may pave the way for a deeper understanding of neuronal synchronization phenomena.

Variations of the phase interaction function

Kuramoto approximated the phase interaction between any two oscillators by its first Fourier component, namely [math]\displaystyle{ \Gamma(\phi) = \sin(\phi) }[/math], where [math]\displaystyle{ \phi = \theta_j - \theta_i }[/math]. Better approximations can be obtained by including higher-order Fourier components,

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

where parameters [math]\displaystyle{ a_i }[/math] and [math]\displaystyle{ b_i }[/math] must be estimated. For example, synchronization among a network of weakly-coupled Hodgkin–Huxley neurons can be replicated using coupled oscillators that retain the first four Fourier components of the interaction function. The introduction of higher-order phase interaction terms can also induce interesting dynamical phenomena such as partially synchronized states, heteroclinic cycles, and chaotic dynamics.

Availability

• pyclustering library includes a Python and C++ implementation of the Kuramoto model and its modifications. Also the library consists of oscillatory networks (for cluster analysis, pattern recognition, graph coloring, image segmentation) that are based on the Kuramoto model and phase oscillator.

See also

• Master stability function
• Oscillatory neural network
• Phase-locked loop
• Swarmalators