Complexity and Networks: The Kuramoto model - Part I by Konstantinos Efstathiou is licensed under a Creative Commons Attribution-ShareAlike 3.0 Unported License.

Collective synchronization is one of the most common phenomena in nature. Some examples include

- flashing fireflies,
- circadian rhythm induced by the synchronization of pacemaker cells in the suprachiasmatic nucleus,

Arthur Winfree proposed in 1967 to study synchronization through the study of interacting limit-cycle oscillators.

Yoshiki Kuramoto in 1975 provided a heuristic analysis of a simple coupled oscillator network that has been since then the basis for further research. This will be the topic of today’s lecture.

Suggested reading:

- Arthur Winfree, The Geometry of Biological Time
- Yoshiki Kuramoto, Chemical Oscillations, Waves, and Turbulence

In this, and part of the next, lecture we follow the exposition in:

The starting point is a system of coupled limit-cycle oscillators.

After parameterization of the motion along the limit cycle by a variable \(\theta_i \in \mathbb S^1 = \mathbb R / 2 \pi \mathbb Z\), such that \(\dot\theta_i = \omega_i\) in the case of no coupling, we can write the equations of motion in the following form

\[ \dot \theta_i = \omega_i + \sum_{j=1}^N F_{ij}(\theta_j - \theta_i), \]

for some functions \(F_{ij}\).

The dynamical system above is very difficult to study at such generality. But we can say some things if we make specific choices for \(F_{ij}\).

Kuramoto’s simplification was to assume that

\[ F_{ij}(\theta_j - \theta_i) = \frac{K}{N} \sin(\theta_j - \theta_i). \]

That is, the coupling between oscillators is exactly the same and it has a very simple form.

The parameter \(K \ge 0\) is the **coupling strength**.

Then the equations of motion become

\[ \dot \theta_i = \omega_i + \frac{K}{N} \sum_{j=1}^N \sin(\theta_j - \theta_i). \]

We introduce the complex-valued **order parameter** by

\[ r e^{i \psi} = \frac{1}{N} \sum_{j=1} e^{i \theta_j}, \]

where \(0 \le r \le 1\) is the modulus of the average of the oscillator phases.

**Remark.** It is worth considering two extreme cases. If all the phases are the same at each moment, that is, \(\theta_j(t) = \theta(t)\) for all oscillators, then \[
r e^{i \psi} = e^{i \theta},
\] implying \(r = 1\). If the phases are uniformly distributed over the circle then \[
\sum_{j=1} e^{i \theta_j} = 0,
\] implying \(r = 0\).

Therefore, \(r\) measures how synchronized the oscillators are. At \(r=0\) we have an **incoherent state** (total lack of synchrony), while at \(r=1\) we have a perfectly synchronous state.

Kuramoto analyzed the model without the help of numerical simulations. Having today much easier access to computers makes it much easier to get a grip on the dynamics of the Kuramoto model through numerical simulations, focusing on the behavior of \(r(t)\).

It turns out that there is a critical value \(K_c\) such that:

When \(K < K_c\) the phases become (almost) uniformly distributed over the circle. Then \(r(t)\) fluctuates near \(0\). This indicates that the incoherent state is a global attractor.

When \(K > K_c\) the incoherent state becomes unstable. Then \(r(t)\) increases until it reaches a value around which it fluctuates. The simulations indicate that this value depends only on \(K\).

Using the order parameter we can write

\[ r e^{i (\psi-\theta_i)} = \frac{1}{N} \sum_{j=1} e^{i (\theta_j-\theta_i)}, \]

and consequently, by taking imaginary parts,

\[ r \sin(\psi-\theta_i) = \frac{1}{N} \sum_{j=1} \sin(\theta_j-\theta_i). \]

Therefore,

\[ \dot \theta_i = \omega_i + \frac{K}{N} \sum_{j=1}^N \sin(\theta_j - \theta_i) = \omega_i + K r \sin(\psi - \theta_i). \]

The equation \[ \dot \theta_i = \omega_i + K r \sin(\psi - \theta_i). \]

shows that the time evolution of each oscillator is determined through the interaction with the whole network through the order parameter, representing an interaction with the **mean field**.

We look for **stationary solutions**, that is, solutions with \(\dot r = 0\) and \(\dot \psi = \Omega\), where we further assume that \(\Omega\) is the mean natural frequency of the oscillators.

Passing to a coordinate system that rotates at the same angular velocity \(\Omega\) as \(\psi\) (using the time-dependent coordinate transformation \(\theta_i \mapsto \theta_i + \Omega t\)) \(\psi\) becomes constant; we choose \(\psi = 0\). Then the equations become

\[\begin{equation}\label{eq:mfeq} \dot \theta_i = \omega_i - \Omega - K r \sin \theta_i. \end{equation}\]We replace \(\omega_i - \Omega\) by \(\omega_i\) with the understanding that the mean of the shifted natural frequencies is \(0\).

We will now analyze the ODE \eqref{eq:mfeq}.

If \(|\omega_i| \le K r\) the ODE has two equilibria, one stable and one unstable, obtained by solving \(\dot\theta_i = 0\). This means that such oscillator will very fast reach the stable equilibrium. In the original coordinates this implies that the oscillator rotates with angular velocity \(\Omega\). We call the corresponding oscillator **locked**.

The stable fixed point is given by solving \[ \sin \theta_i^* = \frac{\omega_i}{K r}. \]

and given that the linearized system is

\[ \dot u_i = - K r \cos \theta_i^* u_i, \]

we have stability for the solution \(\theta_i^*\) that satisfies \(|\theta_i^*| \le \frac12 \pi\).

If \(|\omega_i| > K r\) then the oscillator is **drifting**. We denote by \(\rho(\theta, \omega)\) the distribution of the phases of oscillators with the same natural frequency \(\omega\).

Kuramoto’s assumption was that the distribution of drifting oscillators is **stationary**. For this to happen we must have

\[ \rho(\theta,\omega) = \frac{C}{|\dot \theta|} = \frac{C}{|\omega - K r \sin\theta|}. \]

Normalization gives for each \(\omega\) \[ \int_{-\pi}^\pi \rho(\theta,\omega) d\theta = 1, \] and then \[ C = \frac{1}{2\pi} \sqrt{\omega^2 - (Kr)^2}. \]

Here we have denoted by \(g(\omega)\) the distribution of natural frequencies, that is, the probability density that an oscillator has natural frequency \(\omega\). We assume a symmetric unimodal frequency distribution \(g(\omega)\) with mean value \(0\) and such that \(g(-\omega) = g(\omega)\).

First, we note that \[ \int_{-\pi}^{\pi} \int_{|\omega|>Kr} e^{i\theta} \rho(\theta,\omega) g(\omega) d\omega = 0. \] This is a consequence of the symmetry \(g(\omega)=g(-\omega)\) and the symmetry \(\rho(\theta + \pi, −\omega) = \rho(\theta, \omega)\).

Therefore, \[\begin{align*} r e^{i\psi} = r = \int_{-Kr}^{Kr} e^{i\theta^*(\omega)} g(\omega) d\omega, \end{align*}\]where \[ \omega = Kr \, \sin\theta^*. \]

To compute this integral make the transformation \(u = \theta^*(\omega)\), that is, \(\omega = Kr\,\sin u\). Then \[d\omega = K r \cos u \, du\] and

\[\begin{align*} r & = \int_{-Kr}^{Kr} e^{i\theta^*(\omega)} g(\omega) d\omega \\ & = K r \int_{-\pi/2}^{\pi/2} g(K r \sin u) e^{iu} \cos u \, du \\ & = \frac12 K r \int_{-\pi}^{\pi} g(K r \sin u) e^{iu} \cos u \, du. \end{align*}\]To compute the last integral we can use a standard technique from Complex Analysis that allows us to write trigonometric integrals over \([-\pi,\pi]\) as integrals of complex functions over the unit circle. We then have

\[\begin{align*} r & = \frac12 K r \int_{|z|=1} g \left( \frac{Kr}{2i}\left(z-\frac1z\right) \right) z \frac12 \left(z+\frac1z\right) \frac{dz}{iz} \\ & = \frac{1}{4i} \int_{|z|=1} G(z) \frac{z^2+1}{z} dz, \end{align*}\]where \[ G(z) = K r \, g \left( \frac{Kr}{2i}\left(z-\frac1z\right) \right). \]

This is essentially a self-consistency equation involving \(r\).

Note that \(r = 0\), corresponding to the incoherent state, is always a solution. (For \(r=0\) the integral above equals \(-i/(2\gamma^2)\).)

Solving the above equation for \(r\) we find \[\begin{align*} r = \sqrt{1-\frac{2\gamma}{K}} = \sqrt{1-\frac{K_c}{K}}, \end{align*}\]

where we defined \(K_c = 2\gamma\).

The implication of our results is that for \(K < K_c\) the only possibility is that \(r = 0\), while for \(K > K_c\) there are two possibilities. Either \(r = 0\) or \(r = \sqrt{1-K_c/K}\), as shown in the picture below.

Since the numerical simulations for \(K > K_c\) identify only the branch with \(r = \sqrt{1-K_c/K}\) we expect that this branch is stable, while the branch \(r = 0\) is unstable.

The idea here is to consider the limit \(N \to \infty\) and express this by using distributions of oscillators instead of discrete oscillators. We denote by \[ \rho(\theta,t,\omega) \] the probability density that an oscillator of natural frequency \(\omega\) has phase \(\theta\) at time \(t\).

The normalization condition is that \[ \int_0^{2\pi} \rho(\theta,t,\omega) \, d\theta = 1. \]

Since the number of oscillators is conserved we can write a **continuity equation** which can be derived in the same way as the continuity equation for 1-dimensional fluids. The continuity equation is \[
\frac{\partial \rho}{\partial t} + \frac{\partial }{\partial \theta} (\rho \upsilon) = 0,
\] where \(\upsilon\) is the velocity of an oscillator with natural frequency \(\omega\) and phase \(\theta\). Therefore, \[
\upsilon = \dot\theta = \omega + K r \sin (\psi - \theta).
\]

The order parameter can be written in terms of these distributions as \[ r e^{i \psi} = \int_{0}^{2\pi} \int_{-\infty}^\infty \rho(\theta',t,\omega') g(\omega') e^{i\theta'} \,d\omega'\,d\theta'. \]

Therefore, \[ r e^{i (\psi-\theta)} = \int_{0}^{2\pi} \int_{-\infty}^\infty \rho(\theta',t,\omega') g(\omega') e^{i(\theta'-\theta)} \,d\omega'\,d\theta', \] and we finally obtain by considering the imaginary parts that \[ \upsilon = \omega + K \int_{0}^{2\pi} \int_{-\infty}^\infty \rho(\theta',t,\omega') g(\omega') \sin(\theta'-\theta) \,d\omega'\,d\theta'. \]

The full equation for the continuum limit of the Kuramoto model can then be written as \[ \frac{\partial \rho}{\partial t} + \frac{\partial }{\partial \theta} \left[\rho \left( \omega + K \int_{0}^{2\pi} \int_{-\infty}^\infty \rho(\theta',t,\omega') g(\omega') \sin(\theta'-\theta) \,d\omega'\,d\theta' \right)\right] = 0. \]

Stationary distributions do not change in time. Therefore, \(\partial\rho/\partial t=0\) and \[ \frac{\partial }{\partial \theta} (\rho \upsilon) = 0, \] implying \[ \rho \upsilon = C(\omega). \] If \(C(\omega) \ne 0\), then \(\rho = C(\omega)/\upsilon\). If \(C(\omega) = 0\) then \(\rho\) is a delta-function at \(\upsilon=0\).

A stationary distribution is the uniform distribution \[ \rho(\theta,\omega) = \frac{1}{2\pi}, \] which can be easily seen to satisfy the continuum limit Kuramoto equation and which corresponds to the incoherent state.

Consider the perturbation of the incoherent state given by \[ \rho(\theta,t,\omega) = \frac{1}{2\pi} + \varepsilon \, \eta(\theta,t,\omega). \]

We write a Fourier series for \(\eta\) as \[ \eta(\theta,t,\omega) = \sum_{n=1}^\infty c_n(t,\omega) e^{in\theta} + \bar{c}_n(t,\omega) e^{-in\theta}. \]

Then \[ \int_{0}^{2\pi} \eta(\theta',t,\omega') \sin(\theta'-\theta) \,d\theta' = -2 \pi \mathrm{Im}(c_1 e^{i\theta}). \]

Subsequently, we find that \[ \frac{\partial c_1}{\partial t} + i \omega c_1 - \frac{K}{2} \int_{-\infty}^\infty c_1(t,\omega') g(\omega') d\omega' = 0, \] and, for \(n \ge 2\), \[ \frac{\partial c_n}{\partial t} + i n \omega c_n = 0. \]

Therefore, for \(n \ge 2\) we immediately obtain \[ c_n(t,\omega) = c_n(0,\omega) e^{i n \omega t}. \]

For \(n=1\) we consider the operator \(A\) given by \[ A c := - i \omega c + \frac{K}{2} \int_{-\infty}^\infty c(t,\omega') g(\omega') d\omega', \] so that we can write \[ \frac{\partial c_1}{\partial t} = A c_1. \]

It turns out that the spectrum of \(A\) has a continuous part given by \[ \{ i \omega : \omega \in \mathrm{supp}(g) \}. \] Such eigenvalues imply neutral stability. If we only have eigenvalues on the imaginary axis then the linear stability analysis is not sufficient to determine the stability of the solution.

We look for discrete eigenvalues of \(A\) through the ansatz \[ c(t,\omega) = b(\omega) e^{\lambda t}. \]

Then \[ \lambda b = - i \omega b + \frac{K}{2} \int_{-\infty}^\infty b(\omega') g(\omega') d\omega'. \] Let \[ B = \int_{-\infty}^\infty b(\omega') g(\omega') d\omega'. \] Then \[ b(\omega) = \frac{B}{\lambda + i\omega}. \]

Therefore, we can write \[ \lambda \frac{B}{\lambda + i\omega} = - i \omega \frac{B}{\lambda + i\omega} + \frac{K}{2} \int_{-\infty}^\infty \frac{B}{\lambda+i\omega'} g(\omega') d\omega', \] implying \[ 1 = \frac{K}{2} \int_{-\infty}^\infty \frac{g(\omega')}{\lambda+i\omega'} d\omega'. \] We consider again a Lorentzian distribution \[ g(\omega) = \frac{\gamma}{\pi(\omega^2+\gamma^2)}. \] Therefore, \[ 1 = \frac{K}{2} \int_{-\infty}^\infty \frac{1}{\lambda+i\omega} \frac{\gamma}{\pi(\omega^2+\gamma^2)} d\omega = \frac{K\gamma}{2\pi i} \int_{-\infty}^\infty \frac{1}{\omega-i\lambda} \frac{1}{\omega^2+\gamma^2} d\omega. \]

Using residues we find (closing the integration path with a half-circle on the upper-half complex plane) that for \(\lambda < 0\) we have the contradiction \[ 1 = K \gamma \mathrm{Res}(i\gamma) = \frac{K}{2(\lambda-\gamma)} < 0. \] For \(\lambda > 0\) we find (closing the integration path with a half-circle on the lower-half complex plane) \[ 1 = - K \gamma \mathrm{Res}(-i\gamma) = \frac{K}{2(\lambda+\gamma)}, \] giving \[ \lambda = \frac{K}{2} - \gamma = \frac{K - K_c}{2}. \] This shows that for \(K > K_c\) the incoherent state becomes unstable and small perturbations grow exponentially as \[ \sim \exp\left(\frac{K-K_c}{2} \, t\right). \]