# Kuramoto model

The Kuramoto model (or Kuramoto–Daido model), first proposed by Yoshiki Kuramoto (蔵本 由紀, Kuramoto Yoshiki),[1][2] is a mathematical model used in describing synchronization. More specifically, it is a model for the behavior of a large set of coupled oscillators.[3][4] Its formulation was motivated by the behavior of systems of chemical and biological oscillators, and it has found widespread applications in areas such as neuroscience[5][6][7][8] and oscillating flame dynamics.[9][10] Kuramoto was quite surprised when the behavior of some physical systems, namely coupled arrays of Josephson junctions, followed his model.[11]

The model makes several assumptions, including that there is weak coupling, that the oscillators are identical or nearly identical, and that interactions depend sinusoidally on the phase difference between each pair of objects.

## Definition

In the most popular version of the Kuramoto model, each of the oscillators is considered to have its own intrinsic natural frequency ${\displaystyle \omega _{i}}$, and each is coupled equally to all other oscillators. Surprisingly, this fully nonlinear model can be solved exactly in the limit of infinite oscillators, N→ ∞;[5] alternatively, using self-consistency arguments one may obtain steady-state solutions of the order parameter.[3] The most popular form of the model has the following governing equations:

${\displaystyle {\frac {d\theta _{i}}{dt}}=\omega _{i}+{\frac {1}{N}}\sum _{j=1}^{N}K_{ij}\sin(\theta _{j}-\theta _{i}),\qquad i=1\ldots N}$,

where the system is composed of N limit-cycle oscillators, with phases ${\displaystyle \theta _{i}}$ and coupling constant K.

Noise can be added to the system. In that case, the original equation is altered to

${\displaystyle {\frac {d\theta _{i}}{dt}}=\omega _{i}+\zeta _{i}+{\dfrac {K}{N}}\sum _{j=1}^{N}\sin(\theta _{j}-\theta _{i})}$,

where ${\displaystyle \zeta _{i}}$ is the fluctuation and a function of time. If the noise is considered to be white noise, then

${\displaystyle \langle \zeta _{i}(t)\rangle =0}$ ,
${\displaystyle \langle \zeta _{i}(t)\zeta _{j}(t')\rangle =2D\delta _{ij}\delta (t-t')}$

with ${\displaystyle D}$ denoting the strength of noise.

## Transformation

The transformation that allows this model to be solved exactly (at least in the N → ∞ limit) is as follows:

Define the "order" parameters r and ψ as

${\displaystyle re^{i\psi }={\frac {1}{N}}\sum _{j=1}^{N}e^{i\theta _{j}}}$.

Here r represents the phase-coherence of the population of oscillators and ψ indicates the average phase. Substituting in the equation gives

${\displaystyle {\frac {d\theta _{i}}{dt}}=\omega _{i}+Kr\sin(\psi -\theta _{i})}$.

Thus the oscillators' equations are no longer explicitly coupled; instead the order parameters govern the behavior. A further transformation is usually done, to a rotating frame in which the statistical average of phases over all oscillators is zero (i.e. ${\displaystyle \psi =0}$). Finally, the governing equation becomes

${\displaystyle {\frac {d\theta _{i}}{dt}}=\omega _{i}-Kr\sin(\theta _{i})}$.

## 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 ${\displaystyle \rho (\theta ,\omega ,t)}$. Normalization requires that

${\displaystyle \int _{-\pi }^{\pi }\rho (\theta ,\omega ,t)\,d\theta =1.}$

The continuity equation for oscillator density will be

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

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

${\displaystyle {\frac {\partial \rho }{\partial t}}+{\frac {\partial }{\partial \theta }}[\rho \omega +\rho Kr\sin(\psi -\theta )]=0.}$

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

${\displaystyle re^{i\psi }=\int _{-\pi }^{\pi }e^{i\theta }\int _{-\infty }^{\infty }\rho (\theta ,\omega ,t)g(\omega )\,d\omega \,d\theta .}$

## Solutions for the large N limit

The incoherent state with all oscillators drifting randomly corresponds to the solution ${\displaystyle \rho =1/(2\pi )}$. In that case ${\displaystyle r=0}$, 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

${\displaystyle \rho =\delta \left(\theta -\psi -\arcsin \left({\frac {\omega }{Kr}}\right)\right)}$

for locked oscillators, and

${\displaystyle \rho ={\frac {\rm {normalization\;constant}}{(\omega -Kr\sin(\theta -\psi ))}}}$

for drifting oscillators. The cutoff occurs when ${\displaystyle |\omega |.

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

## 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 ${\displaystyle \omega _{1}=-\omega _{2}}$, and so the system is described exactly by the angle between the two oscillators: ${\displaystyle \Delta \theta =\theta _{1}-\theta _{2}}$. When ${\displaystyle K, the angle cycles around the circle (that is, the fast oscillator keeps lapping around the slow oscillator). When ${\displaystyle K>K_{c}}$, 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 ${\displaystyle \omega _{1},\omega _{2},\omega _{3},K}$, the system has a strange attractor.[13]

## Connection to Hamiltonian systems

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

${\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})}$

After a canonical transformation to action-angle variables with actions ${\displaystyle I_{i}=\left(q_{i}^{2}+p_{i}^{2}\right)/2}$ and angles (phases) ${\displaystyle \phi _{i}=\mathrm {arctan} \left(q_{i}/p_{i}\right)}$, exact Kuramoto dynamics emerges on invariant manifolds of constant ${\displaystyle I_{i}\equiv I}$. With the transformed Hamiltonian

${\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}),}$

Hamilton's equation of motion become

${\displaystyle {\frac {dI_{i}}{dt}}=-{\frac {\partial {\mathcal {H}}'}{\partial \phi _{i}}}=-{\frac {2K}{N}}\sum _{k=1}^{N}{\sqrt {I_{k}I_{i}}}(I_{k}-I_{i})\cos(\phi _{k}-\phi _{i})}$

and

${\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].}$

So the manifold with ${\displaystyle I_{j}=I}$ is invariant because ${\displaystyle {\frac {dI_{i}}{dt}}=0}$ and the phase dynamics ${\displaystyle {\frac {d\phi _{i}}{dt}}}$ becomes the dynamics of the Kuramoto model (with the same coupling constants for ${\displaystyle I=1/2}$). The class of Hamiltonian systems characterizes certain quantum-classical systems including Bose–Einstein condensates.

## Variations of the models

There are a number of types of variations that can be applied to the original model presented above. Some models change to 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[15] (see Transformation and Large N limit above for more info). Network topologies such as rings and coupled populations support chimera states.[16] 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.[17] 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.[18] These are low-amplitude spatially-periodic deviations that emerge from the uniform state (or the wave state) via a Hopf bifurcation.[19] The existence of ripple solutions was predicted (but not observed) by Wiley, Strogatz and Girvan,[20] who called them multi-twisted q-states.

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

A graph with the minimal degree at least ${\displaystyle d_{min}\geq 0.5\ n}$ 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 ${\displaystyle \mu _{c}}$ such that any graph on ${\displaystyle n}$ nodes with minimum degree ${\displaystyle d_{min}\geq \mu _{c}(n-1)}$ must globally synchronise.for ${\displaystyle n}$ large enough. The minimum[20][22] maximum[23] are known to lie between ${\displaystyle 0.6875\leq \mu _{c}\leq 0.75}$.

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

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

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:

${\displaystyle {\frac {d\theta _{i}}{dt}}=\omega _{i}+\sum _{j=1}^{N}a_{ij}\sin(\theta _{j}-\theta _{i}),\qquad i=1\ldots N}$

where ${\displaystyle a_{ij}}$ is a nonzero positive real number if oscillator ${\displaystyle j}$ is connected to oscillator ${\displaystyle i}$. Such model allows for a more realistic study of, e.g., flocking, schooling, and vehicle coordination.[28] 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.[29] Since the Kuramoto model seems to play a key role in assessing synchronization phenomena in the brain,[30] 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 ${\displaystyle \Gamma (\phi )=\sin(\phi )}$, where ${\displaystyle \phi =\theta _{j}-\theta _{i}}$. Better approximations can be obtained by including higher-order Fourier components,

${\displaystyle \Gamma (\phi )=\sin(\phi )+a_{1}\sin(2\phi +b_{1})+...+a_{n}\sin(2n\phi +b_{n})}$,

where parameters ${\displaystyle a_{i}}$ and ${\displaystyle b_{i}}$ 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.[31] The introduction of higher-order phase interaction terms can also induce interesting dynamical phenomena such as partially synchronized states,[32] heteroclinic cycles,[33] and chaotic dynamics.[34]

## 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.

## References

1. ^ Kuramoto, Yoshiki (1975). H. Araki (ed.). Lecture Notes in Physics, International Symposium on Mathematical Problems in Theoretical Physics. Vol. 39. Springer-Verlag, New York. p. 420.
2. ^ Kuramoto Y (1984). Chemical Oscillations, Waves, and Turbulence. New York, NY: Springer-Verlag.
3. ^ a b c Strogatz, Steven H. (2000). "From Kuramoto to Crawford: Exploring the onset of synchronization in populations of coupled oscillators". Physica D. 143 (1–4): 1–20. Bibcode:2000PhyD..143....1S. doi:10.1016/S0167-2789(00)00094-4. S2CID 16668746. Retrieved 2024-01-04.
4. ^ Acebrón, Juan A.; Bonilla, L. L.; Vicente, Pérez; Conrad, J.; Ritort, Félix; Spigler, Renato (2005). "The Kuramoto model: A simple paradigm for synchronization phenomena" (PDF). Reviews of Modern Physics. 77 (1): 137–185. Bibcode:2005RvMP...77..137A. doi:10.1103/RevModPhys.77.137. hdl:2445/12768.
5. ^ a b Bick, Christian; Goodfellow, Marc; Laing, Carlo R.; Martens, Erik A. (2020). "Understanding the dynamics of biological and neural oscillator networks through exact mean-field reductions: a review". Journal of Mathematical Neuroscience. 10 (1): 9. arXiv:1902.05307. doi:10.1186/s13408-020-00086-9. PMC 7253574. PMID 32462281.
6. ^ Cumin, D.; Unsworth, C. P. (2007). "Generalising the Kuromoto model for the study of neuronal synchronisation in the brain". Physica D. 226 (2): 181–196. Bibcode:2007PhyD..226..181C. doi:10.1016/j.physd.2006.12.004. hdl:2292/2666.
7. ^ Breakspear M, Heitmann S, Daffertshofer A (2010). "Generative models of cortical oscillations: Neurobiological implications of the Kuramoto model". Front Hum Neurosci. 4 (190): 190. doi:10.3389/fnhum.2010.00190. PMC 2995481. PMID 21151358.
8. ^ Cabral J, Luckhoo H, Woolrich M, Joensson M, Mohseni H, Baker A, Kringelbach ML, Deco G (2014). "Exploring mechanisms of spontaneous functional connectivity in MEG: How delayed network interactions lead to structured amplitude envelopes of band-pass filtered oscillations". NeuroImage. 90: 423–435. doi:10.1016/j.neuroimage.2013.11.047. hdl:10230/23081. PMID 24321555.
9. ^ Sivashinsky, G.I. (1977). "Diffusional-thermal theory of cellular flames". Combust. Sci. Technol. 15 (3–4): 137–146. doi:10.1080/00102207708946779.
10. ^ Forrester, D.M. (2015). "Arrays of coupled chemical oscillators". Scientific Reports. 5: 16994. arXiv:1606.01556. Bibcode:2015NatSR...516994F. doi:10.1038/srep16994. PMC 4652215. PMID 26582365.
11. ^ Steven Strogatz, Sync: The Emerging Science of Spontaneous Order, Hyperion, 2003.
12. ^ Strogatz, Steven H. (1994). "Norbert Wiener's Brain Waves". In Levin, Simon A. (ed.). Frontiers in Mathematical Biology. Lecture Notes in Biomathematics. Vol. 100. Berlin, Heidelberg: Springer. pp. 122–138. doi:10.1007/978-3-642-50124-1_7. ISBN 978-3-642-50124-1.
13. ^ Maistrenko, Yuri L.; Popovych, Oleksandr V.; Tass, Peter A. (November 2005). "Chaotic Attractor in the Kuramoto Model". International Journal of Bifurcation and Chaos. 15 (11): 3457–3466. Bibcode:2005IJBC...15.3457M. doi:10.1142/S0218127405014155. ISSN 0218-1274.
14. ^ Witthaut, Dirk; Timme, Marc (2014). "Kuramoto Dynamics in Hamiltonian Systems". Phys. Rev. E. 90 (3): 032917. arXiv:1305.1742. Bibcode:2014PhRvE..90c2917W. doi:10.1103/PhysRevE.90.032917. PMID 25314514. S2CID 7510614.
15. ^ Rodrigues, F. A.; Peron, T.K.; Jie, P.; Kurths, J. (2016). "The Kuramoto model in complex networks". Physics Reports. 610 (1): 1–98. arXiv:1511.07139. Bibcode:2016PhR...610....1R. doi:10.1016/j.physrep.2015.10.008. S2CID 119290926.
16. ^ Abrams, D.M.; Strogatz, S.H. (2004). "Chimera states for coupled oscillators". Physical Review Letters. 93 (17): 174102. arXiv:nlin/0407045. Bibcode:2004PhRvL..93q4102A. doi:10.1103/physrevlett.93.174102. PMID 15525081. S2CID 8615112.
17. ^ Kazanci, F.; Ermentrout, B. (2006). "Pattern formation in an array of oscillators with electrical and chemical coupling". SIAM J Appl Math. 67 (2): 512–529. CiteSeerX 10.1.1.140.1020. doi:10.1137/060661041.
18. ^ Heitmann, S.; Gong, P.; Breakspear, M (2012). "A computational role for bistability and traveling waves in motor cortex". Front Comput Neurosci. 6 (67): 67. doi:10.3389/fncom.2012.00067. PMC 3438483. PMID 22973223.
19. ^ Heitmann, S.; Ermentrout, B. (2015). "Synchrony, waves and ripple in spatially coupled Kuramoto oscillators with Mexican hat connectivity". Biological Cybernetics. 109 (3): 1–15. doi:10.1007/s00422-015-0646-6. PMID 25677527. S2CID 18561153.
20. ^ a b Wiley, D.; Strogatz, S.; Girvan, M (2006). "The size of the sync basin". Chaos. 16 (1): 015103. Bibcode:2006Chaos..16a5103W. doi:10.1063/1.2165594. PMID 16599769. S2CID 21173189.
21. ^ Eom, Y.-H.; Boccaletti, S.; Caldarelli, G (2016). "Concurrent enhancement of percolation and synchronization in adaptive networks". Scientific Reports. 7: 015103. arXiv:1511.05468. Bibcode:2016NatSR...627111E. doi:10.1038/srep27111. PMC 4890019. PMID 27251577.
22. ^ Canale, Eduardo A. (2022). "From weighted to unweighted graphs in Synchronizing Graph Theory". arXiv:2209.06362 [math.CO].
23. ^ Kassabov, Martin; Strogatz, Steven H.; Townsend, Alex (2021-07-01). "Sufficiently dense Kuramoto networks are globally synchronizing". Chaos: An Interdisciplinary Journal of Nonlinear Science. 31 (7): 073135. arXiv:2105.11406. Bibcode:2021Chaos..31g3135K. doi:10.1063/5.0057659. ISSN 1054-1500. PMID 34340322. S2CID 235166315.
24. ^ Ling, Shuyang; Xu, Ruitu; Bandeira, Afonso S. (January 2019). "On the Landscape of Synchronization Networks: A Perspective from Nonconvex Optimization". SIAM Journal on Optimization. 29 (3): 1879–1907. arXiv:1809.11083. doi:10.1137/18M1217644. ISSN 1052-6234. S2CID 53067184.
25. ^ Abdalla, Pedro; Bandeira, Afonso S.; Kassabov, Martin; Souza, Victor; Strogatz, Steven H.; Townsend, Alex (2022). "Expander graphs are globally synchronising". arXiv:2210.12788 [math.CO].
26. ^ Sloman, Leila (24 July 2023). "New Proof Shows That 'Expander' Graphs Synchronize". Quanta Magazine.
27. ^ Pantaleone, James (October 2002). "Synchronization of metronomes" (PDF). American Journal of Physics. 70 (10): 992–1000. Bibcode:2002AmJPh..70..992P. doi:10.1119/1.1501118.
28. ^ Dorfler, F.; Bullo, F. (2014). "Synchronization in complex networks of phase oscillators: A survey". Automatica. 50 (6): 1539–1564. doi:10.1016/j.automatica.2014.04.012.
29. ^ Menara, T.; Baggio, G.; Bassett, D.; Pasqualetti, F. (2020). "Stability Conditions for Cluster Synchronizations in Networks of Heterogeneous Kuramoto Oscillators". IEEE Transactions on Control of Network Systems. 7 (1): 302–314. arXiv:1806.06083. doi:10.1109/TCNS.2019.2903914. S2CID 73729229.
30. ^ Cabral, J.; Hugues, E.; Sporns, O.; Deco, G. (2011). "Role of local network oscillations in resting-state functional connectivity". NeuroImage. 57 (1): 130–139. doi:10.1016/j.neuroimage.2011.04.010. PMID 21511044. S2CID 13959959.
31. ^ Hansel, D.; Mato, G.; Meunier, C (1993). "Phase Dynamics for Weakly Coupled Hodgkin-Huxley Neurons". Europhysics Letters. 23 (5): 367–372. Bibcode:1993EL.....23..367H. doi:10.1209/0295-5075/23/5/011. S2CID 250884499.
32. ^ Clusella, Pau; Politi, Antonio; Rosenblum, Michael (2016). "A minimal model of self-consistent partial synchrony". New Journal of Physics. 18 (9): 093037. arXiv:1607.07178. Bibcode:2016NJPh...18i3037C. doi:10.1088/1367-2630/18/9/093037. ISSN 1367-2630. S2CID 126227664.
33. ^ Hansel, D.; Mato, G.; Meunier, C (1993). "Clustering and slow switching in globally coupled phase oscillators" (PDF). Physical Review E. 48 (5): 3470–3477. Bibcode:1993PhRvE..48.3470H. doi:10.1103/physreve.48.3470. PMID 9961005.
34. ^ Bick, C.; Timme, M.; Paulikat, D.; Rathlev, D.; Ashwin, P. (2011). "Chaos in Symmetric Phase Oscillator Networks". Physical Review Letters. 107 (24): 244101. arXiv:1105.2230. Bibcode:2011PhRvL.107x4101B. doi:10.1103/PhysRevLett.107.244101. PMID 22243002. S2CID 16144737.