2.7. The Kuramoto Oscillator

Here, we will introduce the Kuramoto model, a generic phase oscillator model with a wide range of applications [1]. In its simplest form, each Kuramoto oscillator is governed by a non-linear, 1st order ODE:

(1)\[\dot \theta_i = \omega + \sum_j J_{ij} sin(\theta_j - \theta_i),\]

with phase \(\theta\) and intrinsic frequency \(\omega\). The sum represents sinusoidal coupling with all other oscillators in the network with coupling strengths \(J_{ij}\).

In a first step, we’ll consider two coupled Kuramoto oscillators, with an additional extrinsic input \(P(t)\) entering at one of them:

(2)\[\begin{split}\dot \theta_1 &= \omega_1 + P(t) + J_1 sin(\theta_2 - \theta_1), \\ \dot \theta_2 &= \omega_2 + J_2 sin(\theta_1 - \theta_2).\end{split}\]

Below, we will first demonstrate how this model can be used and examined via PyRates.

References

2.7.1. 2 Coupled Kuramoto Oscillators

Here, we use the integrate function imported from PyRates. As a first argument to this function, either a path to a YAML-based model definition or a CircuitTemplate instance can be provided. The function will then compile the model and solve the initial value problem of the above defined differential equations for a time interval from 0 to the given simulation time. This solution will be calculated numerically by a differential equation solver in the backend, starting with a defined step-size. Here, we use the default backend and solver. Furthermore, we provide a step-function extrinsic input to one of the Kuramoto oscillators in a time window from start to stop. This input is defined on a time vector with fixed time steps of size step_size. Check out the arguments of the CircuitTemplate.run() method for a detailed explanation of the arguments that you can use to adjust this numerical procedure.

from pyrates import integrate
import numpy as np

# define simulation time and input start and stop
T = 1.0
step_size = 1e-4
start = 0.2
stop = 0.8

# extrinsic input definition
steps = int(np.round(T/step_size))
I_ext = np.zeros((steps,))
I_ext[int(start/step_size):int(stop/step_size)] = 1.0

# perform simulation
results = integrate("model_templates.oscillators.kuramoto.kmo_2coupled", step_size=step_size, simulation_time=T,
                    outputs={'theta_1': 'p1/phase_op/theta', 'theta_2': 'p2/phase_op/theta'},
                    inputs={'p1/phase_op/s_ext': I_ext}, clear=True)

# plot resulting phases
import matplotlib.pyplot as plt
plt.plot(np.sin(results*2*np.pi))
plt.show()

2.7.2. Kuramoto Order Parameter Dynamics

The Kuramoto order parameter \(z\) of a system of coupled phase oscillators is given by

(3)\[z = \sigma e^{i\phi} = \frac{1}{N} \sum_{j=1}^N e^{i \theta_j},\]

where \(\sigma\) and \(\phi\) are the phase coherence and average phase of the phase oscillators, respectively. In 2008, Ott and Antonsen derived the evolution equations for these order parameters for a system of all-to-all coupled Kuramoto oscillators of the form (1). While the evolution of the average phase is determined by a mere constant, the evolution equation of \(z\) is given by

(4)\[\dot z = z(i\omega - \Delta) - \frac{\bar{J z} z^2 - J z}{2},\]

where \(\bar z\) denotes the complex conjugate of \(z\), and \(\omega\) and \(\Delta\) represent the center and half-width-at-half-maximum of a Lorentzian distribution over the intrinsic frequencies \(\omega_i\) of the individual Kuramoto oscillators (see [2] for a detailed derivation of the mean-field equation). Below, we simulate the dynamics of the dynamics of \(z\) of an all-to-all coupled system of Kuramoto oscillators in response to a step-function input (similar to the simulation above). Note that \(z\) is a complex variable and we plot its absolute value \(|z|\) to receive the coherence of the system over time.

# define simulation time and input start and stop
T = 40.0
step_size = 1e-4
start = 10.0
stop = 30.0

# extrinsic input definition
steps = int(np.round(T/step_size))
I_ext = np.zeros((steps,))
I_ext[int(start/step_size):int(stop/step_size)] = 2.0

# perform simulation
results = integrate("model_templates.oscillators.kuramoto.kmo_mf", step_size=step_size, simulation_time=T,
                    outputs={'z': 'p/kmo_op/z'}, inputs={'p/kmo_op/s_ext': I_ext}, clear=True, solver='scipy',
                    float_precision='complex64')

# plot resulting coherence dynamics
import matplotlib.pyplot as plt
plt.plot(np.abs(results))
plt.show()

As can be seen, the system responded with an increased coherence to the stimulation. Note that we specified float_precision='complex64'. This is important to receice correct results, since it ensures that none of the complex model variables are transformed into real variables, thus losing the imaginary part.

Gallery generated by Sphinx-Gallery