Motivation

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:

• Steven Strogatz, From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators

Equations of motion

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

After parameterization of the motion along the limit cycle by a variable θi∈S1=R/2πZ, such that ˙θi=ωi in the case of no coupling, we can write the equations of motion in the following form

˙θi=ωi+N∑j=1Fij(θj−θi),

for some functions Fij.

The dynamical system above is very difficult to study at such generality. But we can say some things if we make specific choices for Fij.

Kuramoto model

Kuramoto’s simplification was to assume that

Fij(θj−θi)=KNsin(θj−θi).

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

The parameter K≥0 is the coupling strength.

Then the equations of motion become

˙θi=ωi+KNN∑j=1sin(θj−θi).

Order parameter

We introduce the complex-valued order parameter by

reiψ=1NN∑j=1eiθj,

where 0≤r≤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, θj(t)=θ(t) for all oscillators, then reiψ=eiθ, implying r=1. If the phases are uniformly distributed over the circle then ∑j=1eiθ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.

Simulations

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 Kc such that:

1. When K<Kc 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.

2. When K>Kc 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.

Mean field

Using the order parameter we can write

rei(ψ−θi)=1NN∑j=1ei(θj−θi),

and consequently, by taking imaginary parts,

rsin(ψ−θi)=1NN∑j=1sin(θj−θi).

Therefore,

˙θi=ωi+KNN∑j=1sin(θj−θi)=ωi+Krsin(ψ−θi).

The equation ˙θi=ωi+Krsin(ψ−θ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.

Stationary solutions

We look for stationary solutions, that is, solutions with ˙r=0 and ˙ψ=Ω, where we further assume that Ω is the mean natural frequency of the oscillators.

Passing to a coordinate system that rotates at the same angular velocity Ω as ψ (using the time-dependent coordinate transformation θi↦θi+Ωt) ψ becomes constant; we choose ψ=0. Then the equations become

We replace ωi−Ω by ωi with the understanding that the mean of the shifted natural frequencies is 0.

We will now analyze the ODE (1).

If |ωi|≤Kr the ODE has two equilibria, one stable and one unstable, obtained by solving ˙θ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 Ω. We call the corresponding oscillator locked.

The stable fixed point is given by solving sinθ∗i=ωiKr.

and given that the linearized system is

˙ui=−Krcosθ∗iui,

we have stability for the solution θ∗i that satisfies |θ∗i|≤12π.

If |ωi|>Kr then the oscillator is drifting. We denote by ρ(θ,ω) the distribution of the phases of oscillators with the same natural frequency ω.

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

ρ(θ,ω)=C|˙θ|=C|ω−Krsinθ|.

Normalization gives for each ω ∫π−πρ(θ,ω)dθ=1, and then C=12π√ω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 reiψ=1NN∑i=1eiθi=1NNL∑i=1eiθ∗i+1NND∑i=1eiθi.

The next step is to pass to the continuous limit and replace these sums by integrals. We then have reiψ=∫Kr−Kreiθ∗(ω)g(ω)dω+∫π−π∫|ω|>Kreiθρ(θ,ω)g(ω)dωdθ.

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

First, we note that ∫π−π∫|ω|>Kreiθρ(θ,ω)g(ω)dωdθ=0. This is a consequence of the symmetry g(ω)=g(−ω) and the symmetry ρ(θ+π,−ω)=ρ(θ,ω).

Therefore, reiψ=r=∫Kr−Kreiθ∗(ω)g(ω)dω, where ω=Krsinθ∗.

To compute this integral make the transformation u=θ∗(ω), that is, ω=Krsinu. Then dω=Krcosudu and

To compute the last integral we can use a standard technique from Complex Analysis that allows us to write trigonometric integrals over [−π,π] as integrals of complex functions over the unit circle. We then have

where G(z)=Krg(Kr2i(z−1z)).

The last integral can be explicitly computed for a Lorentzian frequency distribution g(ω)=γπ(γ2+ω2). We find G(z)=−4γKrz2πK2r2(z2−1)2−4πγ2z2, therefore, r=iγKr∫|z|=1z(z2+1)πK2r2(z2−1)2−4πγ2z2dz.

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γ2).)

Exercise Prove, using residues, that for r≠0 we have, r=√1+(γKr)2−γKr.

Solving the above equation for r we find r=√1−2γK=√1−KcK, where we defined Kc=2γ.

The implication of our results is that for K<Kc the only possibility is that r=0, while for K>Kc there are two possibilities. Either r=0 or r=√1−Kc/K, as shown in the picture below.

Since the numerical simulations for K>Kc identify only the branch with r=√1−Kc/K we expect that this branch is stable, while the branch r=0 is unstable.