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

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

2.6.1. References

1
  1. Kuramoto (1991) Collective synchronization of pulse-coupled oscillators and excitable units. Physia D: Nonlinear Phenomena 50(1): 15-30.

2

Ott and Antonsen (2008) Low dimensional behavior of large systems of globally coupled oscillators. Chaos 18(3): 1054-1500.

2.6.1.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.6.1.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')

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

As can be seen, the system engaged in synchronized oscillations within the input period.

Total running time of the script: ( 0 minutes 0.000 seconds)

Gallery generated by Sphinx-Gallery