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}^N 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}^N e^{i (\theta_j-\theta_i)}, \]

and consequently, by taking imaginary parts,

\[ r \sin(\psi-\theta_i) = \frac{1}{N} \sum_{j=1}^N \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}. \]

We can write the expression giving the order parameter as two sums, one involving locked oscillators, and one involving drifting oscillators. We have \[\begin{align*} r e^{i\psi} & = \frac{1}{N} \sum_{i=1}^N e^{i\theta_i} \\ &= \frac{1}{N} \sum_{i=1}^{N_L} e^{i\theta_i^*} + \frac{1}{N} \sum_{i=1}^{N_D} e^{i\theta_i}. \end{align*}\]

The next step is to pass to the continuous limit and replace these sums by integrals. We then have \[\begin{align*} r e^{i\psi} = \int_{-Kr}^{Kr} e^{i\theta^*(\omega)} g(\omega) d\omega + \int_{-\pi}^{\pi} \int_{|\omega|>Kr} e^{i\theta} \rho(\theta,\omega) g(\omega) \, d\omega \, d\theta. \end{align*}\]

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 \, d\theta = 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). \]

The last integral can be explicitly computed for a Lorentzian frequency distribution \[\begin{align*} g(\omega) = \frac{\gamma}{\pi(\gamma^2+\omega^2)}. \end{align*}\] We find \[\begin{align*} G(z) = -\frac{4 \gamma K r z^2}{\pi K^2 r^2 (z^2-1)^2 - 4\pi \gamma^2 z^2}, \end{align*}\] therefore, \[\begin{align*} r & = i \gamma K r \int_{|z|=1} \frac{z (z^2+1)}{\pi K^2 r^2 (z^2-1)^2 - 4\pi \gamma^2 z^2} dz. \end{align*}\]

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)\).)

**Exercise** Prove, using residues, that for \(r \ne 0\) we have, \[\begin{align*}
r
= \sqrt{1+\left(\frac{\gamma}{Kr}\right)^2}-\frac{\gamma}{Kr}.
\end{align*}\]

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.