Simulations of Kuramoto Model

Apr 2019
I'm trying to understand and computationally model the Kuramoto Model. This is a mathematical model used to describe synchronization (phase locking of a set of $N$ oscillators). The governing equations: $$\frac{d \theta_i}{dt} = \omega_i + \frac{K}{N}\sum_{j=1}^{N}\sin(\theta_j - \theta_i)~~~~~~~~i=1...N$$
where the system is composed of $N$ oscillators, with phases $\theta_i$ and coupling $K$.

As I understand, these phases $\theta_i$ refer to the respective phases of each oscillator. For an individual harmonic oscillator this would be $\phi$ in $$x(t) = A \cos(\omega t + \phi).$$

Consider this animation which shows various levels of phase locking by plotting the phases of the respective oscillators as a function of time.

Question: I'm trying to find a viable numerical scheme used to produce this type of animation. Would it be viable to use the transformation $$re^{i \psi} = \frac{1}{N} \sum_{j=1}^{N}e^{\theta_j}$$ which allows a reformulation of the governing equations to $$\frac{d \theta_i}{d t} = \omega_i - Kr \sin(\theta_i).$$ As described in the wiki entry Kuramoto model. Then using some initial phase values $\theta_i$, together with the reformulated governing equations, we could map out the evolution of each individual phase as in the animation. Is this a viable numeric scheme for producing the animation which shows phase locking?

Thanks for your assistance.