Stochastic Kuramoto oscillators with discrete phase states

TL;DR

This paper introduces a stochastic generalization of the Kuramoto model with discrete phase states, analyzing how phase discretization affects synchronization and oscillation quality through analytical and numerical methods.

Contribution

It presents a novel stochastic Kuramoto model with discrete phases, bridging continuous and discrete oscillatory systems, and explores its synchronization properties.

Findings

Key observables show extrema at certain discretization levels.

Steady-state synchrony varies with phase discretization.

Model converges to classical Kuramoto in the continuous limit.

Abstract

We present a generalization of the Kuramoto phase oscillator model in which phases advance in discrete phase increments through Poisson processes, rendering both intrinsic oscillations and coupling inherently stochastic. We study the effects of phase discretization on the synchronization and precision properties of the coupled system both analytically and numerically. Remarkably, many key observables such as the steady-state synchrony and the quality of oscillations show distinct extrema while converging to the classical Kuramoto model in the limit of a continuous phase. The phase-discretized model provides a general framework for coupled oscillations in a Markov chain setting.

Equations

$$\frac{\mathrm{d}{\varphi }_i}{\mathrm{d}t}={\omega }_i+\kappa \sum _{j=1}^N{c}_{ij}\Gamma ({\varphi }_j-{\varphi }_i),$$

$$\epsilon \frac{\partial P}{\partial t}$$

$$P(\phi ,t|{\phi }^{\prime },0)=\mathrm{Poisson}(\omega t/\epsilon ,\phi -{\phi }^{\prime }),$$

$$\tilde{\omega }=\frac{\sin \epsilon }{\epsilon }\omega ,k=\frac{1-\cos \epsilon }{\epsilon }\omega .$$

$$Q=\frac{1}{2\pi \tan (\epsilon /2)}.$$

$$Q=\frac{m}{2{\pi }^2}-\frac{1}{6m}+\mathcal{O}(m^{-3})$$

$$\frac{\mathrm{d}\varphi }{\mathrm{d}t}=\omega +\sqrt{\frac{2\pi \omega }{m}}\eta (t),$$

$$\epsilon \frac{\partial P}{\partial t}$$

$${\hat{\omega }}_i$$

$${\hat{\phi }}_i^{\pm }P({\phi }_1,\dots ,{\phi }_N,t)=P({\phi }_1,\dots ,{\phi }_i\pm 1,\dots ,{\phi }_N,t).$$

$$\begin{array}{rl}{\hat{\Gamma }}_{ij}(\phi )& =\frac{1+\Gamma (\epsilon [\phi +1])}{2}{\hat{\phi }}_i^{-}-1\\ & +\frac{1-\Gamma (\epsilon [\phi -1])}{2}{\hat{\phi }}_i^{+}.\end{array}$$

$$\displaystyle\phi_{i}(t)=\Phi_{i}(t)+\sqrt{\frac{2\pi}{m}}\xi_{i}(t)+\mathcal{O}(m^{-1})\ ,\bigg{.}$$

$$\frac{\mathrm{d}{\Phi }_i}{\mathrm{d}t}$$

$$\frac{\mathrm{d}{\xi }_i}{\mathrm{d}t}$$

$${\mu }_i=|{\omega }_i|+\kappa c_i,$$

$$\Gamma (\varphi )=\sin (\varphi -{\varphi }_0),$$

$$\Psi (t)=\frac{1}{N}\sum _{i=1}^N{x}_i(t),$$

$$\rho (t)\equiv ⟨|\Psi (t){|}^2⟩=\frac{1+\mathrm{Re}\chi (t)}{2}.$$

$$R$$

$${\tau }_{\mathrm{h}}$$

$${\tau }_{\nu }$$

$$\gamma =\frac{|{\omega }_1-{\omega }_2|}{2\kappa \cos {\varphi }_0}.$$

$$R=\frac{1+\mathrm{sign}(\gamma )\sqrt{1-{\gamma }^2}}{2}.$$

$$R$$

$${\Lambda }_n=-\underset{i=n}{\overset{m-1}{K}}{\lambda }_i\equiv -\frac{1}{{\lambda }_n+\frac{1}{{\lambda }_{n+1}+\frac{1}{\ddots +\frac{1}{{\lambda }_{m-1}}}}}$$

$${\lambda }_n$$

$$\epsilon \frac{\partial P}{\partial t}$$

$$\begin{array}{rl}{\hat{\Gamma }}_i^{\mathrm{MF}}(\phi ,\Psi )& =\sum _{\sigma =\pm }{\left[\mathrm{Im}\left\{\frac{N\Psi {\mathrm{e}}^{-\mathrm{i}\epsilon \phi }-1}{N-1}{\mathrm{e}}^{\mathrm{i}(\sigma \epsilon -{\varphi }_0)}\right\}\right]}_{\sigma }{\hat{\phi }}_i^{-\sigma }\\ & -1,\end{array}$$

$$\frac{\mathrm{d}{\Phi }_i}{\mathrm{d}t}$$

$$\displaystyle\begin{split}\frac{\mathrm{d}\xi_{i}}{\mathrm{d}t}&=\frac{\kappa}{1-N^{-1}}\Big{\{}\tilde{r}\cos(\tilde{\psi}-\Phi_{i}-\phi_{0})\\ &\qquad-r\xi_{i}\cos(\psi-\Phi_{i}-\phi_{0})\Big{\}}+\sqrt{\mu_{i}}\eta_{i}(t)\ ,\end{split}$$

$$\displaystyle r\mathrm{e}^{\mathrm{i}\psi}=\frac{1}{N}\sum_{j}\mathrm{e}^{\mathrm{i}\Phi_{j}}\ ,\qquad\tilde{r}\mathrm{e}^{\mathrm{i}\tilde{\psi}}=\frac{1}{N}\sum_{j}\xi_{j}\mathrm{e}^{\mathrm{i}\Phi_{j}}\ .\Bigg{.}$$

Full text

I Introduction

The dialectic of synchronization has become a powerful conceptual tool in theoretical physics—rooted in the description of coupled oscillators and clocks Strogatz and Stewart (1993), it has been extended to phenomena that bear only structural resemblance to coupled oscillators such as the collective behavior of bird flocks Ha et al. (2010) and magnetic systems Flovik et al. (2016). Hence, it is not surprising that probably the most prominent theoretical paradigm for synchronization, the celebrated Kuramoto model of coupled phase oscillators and its multifarious variants Kuramoto (1984a, b); Acebrón et al. (2005); Rodrigues et al. (2016), have been applied to problems as different as neutrino oscillations Pantaleone (1998), embryonic body axis segmentation Giudicelli et al. (2007); Morelli et al. (2009); Jörg et al. (2016), electric power grids Filatrella et al. (2008); Rohden et al. (2012); Dörfler and Bullo (2012), epileptic seizures Schmidt et al. (2014), and quantum entanglement Witthaut et al. (2017). The Kuramoto model is a time-continuous, phase-continuous system of coupled differential equations Kuramoto (1984b); Acebrón et al. (2005); Rodrigues et al. (2016),

[TABLE]

where $\phi_{i}$ is the phase of oscillator $i=1,\ldots,N$ and $\omega_{i}$ is its intrinsic frequency, $\kappa$ is the coupling strength, $\Gamma$ is a $2\pi$-periodic coupling function, and $c_{ij}$ is the coupling topology matrix where $c_{ij}>0$ indicates that the dynamics of oscillator $i$ couples to the dynamics of oscillator $j$ and $c_{ij}=0$ otherwise. For appropriate choices of $\Gamma$ and $c_{ij}$, the coupling term alters the dynamic frequency $\mathrm{d}\phi_{i}/\mathrm{d}t$ in such a way that the system tends to synchronize, given that coupling can overcome the spread in frequencies Strogatz (2000).

Whether a phase-continuous model is a viable description depends on the system at hand. Biochemical oscillators, for instance, operate through chemical and/or genetic feedbacks between different molecule species and are often characterized by small numbers of molecules which are subject to fluctuations Morelli and Jülicher (2007); Zwicker et al. (2010); Suvak and Demir (2012); Webb et al. (2016); Lengyel and Morelli (2017). Another prominent example from biology is the cell cycle, which, while going through well-defined states, can exhibit considerable period variations Weber et al. (2014). Often, it is desirable to represent such processes on a coarse-grained level, e.g., by Markov chain models, when only their core features are to be retained. This is especially interesting if coupled oscillatory processes are part of a more complex system involving interactions with non-oscillatory parts. The latter is often the case in biology, where periodic processes interact with cell fates, intercellular signaling systems, and/or tissue growth Sugimoto et al. (2012); Brown (2014); Jörg et al. (2016). In recent years, there has been an extensive interest in the behavior of discrete-state models of uncoupled and coupled oscillators Prager et al. (2003); Wood et al. (2006, 2007); Fernandez and Tsimring (2008); Tönjes and Kori (2011); Suvak and Demir (2012); Pinto et al. (2014); Escaff et al. (2014); Barato and Seifert (2016). Recently, for instance, the question has been investigated whether discrete-state models can capture the behavior of the noisy Kuramoto model with all-to-all coupling and homogeneous frequencies Escaff et al. (2016).

In this paper, we study a generalization of the Kuramoto model in which each oscillator transitions between discrete phase states with defined transition rates. This renders all parts of the model inherently stochastic, including the coupling dynamics between oscillators. We investigate the effects of phase discretization on the dynamics of systems with homogeneous and inhomogeneous frequencies, in particular their synchronization behavior, their phase-coherence, and their period fluctuations. In Section II, we introduce the description of a single phase-discretized oscillator and characterize its stochastic properties such as its effective frequency and its quality factor. In Section III, we introduce a stochastic generalization of the coupled Kuramoto model with arbitrary coupling topology and coupling function and discuss variants of this generalization. In Section IV, we study the case of two coupled oscillators and investigate the effects of coupling on synchronization and precision both analytically and numerically. In Section V, we consider the case of many oscillators with homogeneous frequencies and present numerical results on their synchronization behavior and their collective precision. Moreover, we study the onset of synchronization in a system with inhomogeneous frequencies. Finally, in Section VI, we briefly summarize our results, discuss their relevance, and suggest directions for further studies.

II A single phase-discretized oscillator

We start by considering a single phase-discretized oscillator. We discretize the phase interval $[0,2\pi)$ into $m$ states and allow the oscillator to advance by discrete phase increments of size $\varepsilon=2\pi/m$, so that its state is given by the discrete phase variable $\varphi\in\mathds{Z}$ (Fig. 1a). The discrete state $\varphi$ is associated with a phase $\phi=\varepsilon\varphi\in\mathds{R}$ and the corresponding oscillatory signal $x=\exp(\mathrm{i}\phi)$.

The stochastic dynamics of the oscillator is governed by a master equation for the probability $P=P(\varphi,t)$ that the oscillator has the discrete phase $\varphi$ at time $t$. Introducing a transition frequency $\omega\geq 0$, we describe the transition $\varphi\to\varphi+1$ as a Poisson process with transition rate $\omega/\varepsilon$ for a given discretization $m$ (Fig. 1a). This ensures that the average duration of one revolution is given by $2\pi/\omega$. The corresponding master equation is given by

[TABLE]

The solution to Eq. (2) for the initial condition $P(\varphi,0)=\delta_{\varphi\varphi^{\prime}}$ is a Poisson distribution van Kampen (2011),

[TABLE]

where $\mathrm{Poisson}(\lambda,n)=\lambda^{n}\smash{\mathrm{e}^{-\lambda}}\Theta(n)/n!$ with $\Theta$ being the Heaviside function. Fig. 1b shows examples of stochastic trajectories for different phase discretizations $m$, obtained from a standard stochastic simulation algorithm Gillespie (1977).

To determine the dynamical frequency of the oscillator and the frequency fluctuations introduced by stochasticity, we compute the temporal autocorrelation function $G(t)=\langle x(t)x^{*}(0)\rangle$ of the associated oscillatory signal, where the star denotes the complex conjugate. For the system specified by Eq. (2), it assumes the form $G(t)=\exp(\mathrm{i}\tilde{\omega}t-kt)$, where the effective frequency $\tilde{\omega}$ and the decorrelation rate $k$ are given by

[TABLE]

Note that both the effective frequency and the decorrelation rate are proportional to $\omega$, which is the only (inverse) time scale in the system. Notably, the effective frequency $\tilde{\omega}$ is systematically smaller than the transition frequency $\omega$. This difference is due to a ‘stroboscopic’ effect: Starting from a defined state $\varphi^{\prime}$ at time 0, Eq. (3) implies that the discrete phase increment $\Delta\varphi=\varphi-\varphi^{\prime}$ at $t>0$ is Poisson-distributed with mean and variance depending on $m$ through $\varepsilon$. For coarse phase discretizations $m$, the tail of the Poisson distribution can considerably extend into regions with phase increments $\Delta\varphi>m$; that is, for a given elapsed time interval, there is a non-vanishing probability for the oscillator to advance by more than one complete cycle. However, since the oscillatory signal $x$ is $m$-periodic in $\varphi$, a phase increment $\Delta\varphi$ larger than $m$ leads to the same signal as the increment $(\Delta\varphi\operatorname{mod}m)<m$, implying that in this case, the oscillator seemingly goes more slowly. As expected, in the limit of a continuous phase $\varepsilon\to 0$, we recover $\tilde{\omega}\to\omega$ and $k\to 0$.

Frequency fluctuations are commonly characterized by the quality factor $Q=\tilde{\omega}/(2\pi k)$ Pikovsky and Kurths (1997); Morelli and Jülicher (2007), which is independent of the absolute frequency scale and corresponds to the number of oscillations over which the oscillatory signal stays correlated,

[TABLE]

For large phase discretizations $m$, the quality factor quickly approaches the asymptotic behavior

[TABLE]

and becomes effectively linear in $m$. Hence, even to achieve a very low quality factor of $Q=1$, about $m=20$ internal states are required.

To make the connection to the Kuramoto model Eq. (1), we derive the Langevin equation describing the system in the large-$m$ limit by a system size expansion van Kampen (2011), see Appendix A. This yields

[TABLE]

where $\eta$ is Gaussian white noise with $\langle\eta(t)\rangle=0$ and $\langle\eta(t)\eta(t^{\prime})\rangle=\delta(t-t^{\prime})$. The noise strength grows with the transition frequency $\omega$, a reflection of the fact that the quality factor—which in the case of the Langevin equation (7) is given by the linear term in Eq. (6)—only depends on the phase discretization $m$ but not the transition frequency $\omega$, which provides the only time scale in the system and, therefore, cannot affect any dimensionless physical quantity.

III Description of oscillator coupling

While the formalization of an uncoupled phase-discretized oscillator seems straightforward, introducing coupling to such a system opens a plethora of different possibilities, even if coupling processes are constrained to a set of Poisson processes running in parallel. A general formulation of oscillator coupling inevitably introduces backward jumps of the phase since the coupling strength may exceed the intrinsic frequency of an oscillator and may therefore lead to a negative dynamic frequency. While different stochastic formulations can produce the same mean-field dynamics and/or the same phase-continuous limit $m\to\infty$, fluctuations depend on the details of the stochastic dynamics, e.g., whether for each oscillator, (i) forward/backward jumps of the discrete phase are independent processes running in parallel or (ii) whether only forward or only backward processes can occur depending on the phase relation to other oscillators. We give a brief discussion of different possibilities in Appendix B.

III.1 Stochastic Kuramoto model with discretized phases

With these different possibilities in mind, we can now write a specific stochastic formulation of coupled phase-discretized oscillators. The probability $P=P(\varphi_{1},\ldots,\varphi_{N},t)$ of $N$ oscillators with discrete states $\varphi_{1},\ldots,\varphi_{N}\in\mathds{Z}$ is governed by the master equation

[TABLE]

where the operators $\hat{\omega}_{i}$ and $\hat{\Gamma}_{ij}$ describe the stochastic dynamics of intrinsic oscillations and coupling, respectively. As in the Kuramoto model Eq. (1), $\kappa\geq 0$ denotes the coupling strength and $c_{ij}$ is the adjacency matrix determining the coupling topology. The intrinsic frequency operator $\hat{\omega}_{i}$ assumes the generic form

[TABLE]

where $\omega_{i}$ is the transition frequency of oscillator $i$ and where we have used the ladder operator notation $\smash{\hat{\varphi}_{i}^{\pm}}$, defined by van Kampen (2011)

[TABLE]

For a given $2\pi$-periodic coupling function $\Gamma(\phi)$ taking values between $-1$ and $1$, we define the coupling operator as

[TABLE]

This specific formulation of the coupling term corresponds to a biased discrete diffusion process on the discretized phase space, where the bias dynamically depends on the phase difference through the coupling function $\Gamma(\phi)$; i.e., depending on the phase difference, one of the processes $\varphi_{i}\to\varphi_{i}+1$ or $\varphi_{i}\to\varphi_{i}-1$ is favored. Note that the coupling operator given by Eq. (11) does not explicitly depend on the index $j$ of the sending oscillator as compared to other discretization schemes, see Appendix B.

Note that in Eq. (8), the expression in parentheses formally resembles the r.h.s. of the Kuramoto model Eq. (1) with parameters and functions promoted to Liouville operators. In general, the phase discretization leads to two major differences to the classical Kuramoto model Eq. (1): (i) oscillator dynamics is now inherently stochastic and (ii) the coupling function is sampled at discrete readout points determined by the phase discretization.

III.2 Linear noise approximation

To establish a connection to the classical Kuramoto model Eq. (1), we carry out a system size expansion of the system described by Eqs. (8–11), formally interpreting the phase discretization $m$ as the system size. This enables to write a linear noise approximation (LNA) for the corresponding stochastic governing equations for the physical phases $\phi_{i}=\varepsilon\varphi_{i}$. Details on the derivation are given in Appendix A. The resulting linear noise approximation for the physical phase is given by

[TABLE]

where $\Phi_{i}$ is a ‘macroscopic’ phase variable obeying the deterministic Kuramoto dynamics

[TABLE]

cf. Eq. (1), and $\xi_{i}$ is a random variable which is governed by the Langevin equation

[TABLE]

where $\Gamma^{\prime}$ is the derivative of the coupling function, $\eta_{i}$ is Gaussian white noise with $\langle\eta_{i}(t)\rangle=0$ and $\langle\eta_{i}(t)\eta_{j}(t^{\prime})\rangle=\delta_{ij}\delta(t-t^{\prime})$, and $\mu_{i}$ is the effective noise strength for oscillator $i$, given by

[TABLE]

where $c_{i}=\sum_{j}c_{ij}$ is the total coupling weight of the oscillators coupled into oscillator $i$. Eq. (15) illustrates that in our formulation of the phase-discretized system, noise has two sources: the intrinsic oscillatory dynamics, as indicated by the intrinsic transition frequency $\omega_{i}$ and already shown in Eq. (7), but also coupling which, as an inherently stochastic process, inevitably contributes noise to the system. For a given oscillator $i$, coupling to each connected oscillator $j$ is an independent process; therefore, the total contribution from coupling to its noise strength is proportional to the total coupling weight $c_{i}$. Hence, the net effect of coupling on the synchronization and precision properties of the coupled system are not immediately obvious. Note that in the limit of a continuous phase $m\to\infty$, the phases behave as $\phi_{i}\to\Phi_{i}$ and the classical Kuramoto model Eq. (1) for the physical phases $\phi_{i}$ is recovered.

Note also that for the master equation of the single oscillator, Eq. (2), which has state-independent transition rates, the derived linear noise approximation is exact up to second order in the moments Grima (2015). On the other hand, it is not a priori obvious under which circumstances the linear noise approximation Eq. (12–15) is a good approximation for the full model Eqs. (8–11), as it involves (i) nonlinearly state-dependent transition rates and (ii) entails an expansion of the coupling function $\Gamma$, suggesting that its validity is constrained to the vicinity of states for which $\Gamma$ is approximately linear around the occurring phase differences. In Section V, we demonstrate its effectiveness in describing steady-state properties by numerical simulations.

IV Dynamics of two coupled oscillators

To gain some insights into the stochastic behavior of the coupled system, we first investigate the simplest case of two coupled oscillators without self-coupling, $N=2$ and $c_{ij}=1-\delta_{ij}$. For concreteness, we consider the generic class of coupling functions of the Kuramoto–Sakaguchi type in the following Sakaguchi et al. (1988),

[TABLE]

where $\phi_{0}\in[0,2\pi)$ is a constant phase shift. First, we address the transient behavior of the oscillators approaching synchrony, followed by an analysis of the synchronization and precision properties of the steady state.

IV.1 Synchronization transient

The time-dependent phase coherence of the two oscillators can be monitored via the Kuramoto order parameter, defined by Kuramoto (1984b)

[TABLE]

where $N$ is the number of oscillators and $x_{i}(t)=\mathrm{e}^{\mathrm{i}\varepsilon\varphi_{i}(t)}$ is the oscillatory signal associated with oscillator $i$, as before. Usually, one considers the magnitude $|\Psi|$, which takes values from [math] to $1$ with $|\Psi|=1$ indicating perfect phase coherence and $|\Psi|<1$ indicating the existence of phase lags between oscillators. Here we focus on the squared magnitude $|\Psi|^{2}$, which basically has the same interpretation but simpler analytical properties. For two oscillators, the cross correlation $\chi(t)=\langle x_{1}(t)x_{2}^{*}(t)\rangle$ contains the expectation value of $|\Psi|^{2}$ in its real part,

[TABLE]

Since $|\Psi|^{2}$ is bounded, a value of $\rho$ close to 1 indicates not only a small average phase difference but also small fluctuations in the phase difference. Fig. 2 shows both $\rho(t)$ as well as $\operatorname{Im}\chi(t)$ (which together carry the same information as the full cross correlation $\chi$) for different phase discretizations $m$ for two oscillators with unequal frequencies and an initially maximally desynchronized state. After an initial transient, the system approaches a steady state with a constant order parameter and cross correlation

[TABLE]

which depend on $m$.

Even though coarser phase discretizations typically entail a lower degree of synchrony at steady-state, such systems tend to initially synchronize faster than system with a finer discretization (see Fig. 2a). To illustrate this behavior, we define, for a system starting from the completely desynchronized state with maximum phase difference, as two complementary quantities the time $\tau_{\mathrm{h}}$ it takes for the order parameter to reach the absolute value $1/2$ and the time $\tau_{\nu}$ it takes to reach a relative fraction $\nu$ of the steady-state order parameter $R$,

[TABLE]

Fig. 3 shows the synchronization times as a function of the phase discretization $m$ for different values of $\nu$ and reveals an interesting behavior: The time to reach the absolute order parameter $\rho=1/2$ tends to decrease for coarser phase discretizations even though the transition frequencies $\omega_{i}$ and the coupling strength $\kappa$ are kept constant (Fig. 3a). In contrast, the time to reach a relative fraction of the steady-state order parameter attains a distinct maximum for finite discretizations (Fig. 3b). Therefore, coarser phase discretizations can facilitate faster initial synchronization even though they eventually reach a smaller phase-coherence and take a longer time reach the vicinity of their steady state.

IV.2 Steady-state phase coherence

How does the steady-state phase coherence depend on the phase discretization? And how does this compare to the synchronized state of two coupled Kuramoto oscillators with detuning? Let us briefly recapitulate some results from the classical Kuramoto model Adler (1946); Schuster and Wagner (1989). There, the system assumes a phase-locked steady state if coupling is strong enough to overcome the frequency difference between the oscillators, that is, if $|\gamma|<1$ where

[TABLE]

In this case, the order parameter is given by

[TABLE]

Hence, in terms of the intrinsic frequences, the order parameter is determined by the absolute frequency difference $|\omega_{1}-\omega_{2}|$ in a monotonic way. For $|\gamma|>1$, both oscillators phase-drift with respect to each other and the time average of the order parameter is $1/2$.

In the case of the phase-discretized model, nonlinear coupling combines with stochasticity and therefore, an analysis is more involved. Nevertheless, an exact solution for the steady-state order parameter $R$ and the cross correlation $X$, Eqs. (19), can be constructed, see Appendix C for a derivation. Without loss of generality, we consider the case $\omega_{1}\geq 0$, $\omega_{2}\geq 0$ for which the resulting order parameter is given by

[TABLE]

where $\Lambda_{n}$ can be represented as the continued fraction

[TABLE]

with

[TABLE]

Interestingly, the order parameter $R$ depends not only on the frequency difference $\omega_{1}-\omega_{2}$ but also on the frequency sum $\omega_{1}+\omega_{2}$ through $\lambda_{n}$. This reflects the fact that in the stochastic system, the degree of noise depends on the frequency scale (cf. Eq. (7) and the discussion below). Due to its combinatorial complexity, the exact solution given by Eq. (24–26) is somewhat opaque; therefore, we give a few explicit expressions for small phase discretizations $m$ in Appendix C.

Fig. 4 shows the order parameter $R$ and the imaginary part of the cross correlation ${X}$ as a function of the phase discretization $m$ for different frequency detunings and coupling strengths, both from numerical simulations of stochastic trajectories (dots) and the exact solutions given by Eqs. (24) and (43) (solid lines). For many generic parameter combinations, the order parameter monotonically increases with finer phase discretizations. However, in a few cases, the behavior of the order parameter and the cross correlation show some remarkable features. First, even at coupling strengths below the classical critical value that ensures $|\gamma|<1$ we detect partial synchrony, i.e., an order parameter $R>1/2$ (bright curve in Fig. 4c, corresponding to $\gamma=1.25$), indicating that the system spends a larger time in regions with small phase differences. Second, while in all cases the order parameter approaches the Kuramoto value in the limit $m\to\infty$, the convergence is not always monotonic. In fact, there are phase discretizations $m$ for which the degree of partial synchrony becomes maximal. This is exemplified by the bright curve in Fig. 4c and in Fig. 5, where the order parameter is displayed for different coupling strengths and up to very fine phase discretizations. The curves below the critical coupling strength $\kappa_{\mathrm{c}}=|\omega_{1}-\omega_{2}|/2$ exhibit a non-monotonic behavior with a distinct maximum for a finite phase discretization.

This behavior can be illuminated as follows: In the deterministic case $m\to\infty$, the phase difference $\psi=\phi_{1}-\phi_{2}$ of both oscillators is governed by the Adler equation $\mathrm{d}\psi/\mathrm{d}t=-\mathrm{d}v/\mathrm{d}\psi$ with $v(\psi)=-(\omega_{1}-\omega_{2})\psi-2\kappa\cos\psi$ Adler (1946), where for simplicity, we have considered the case of zero coupling phase shift, $\phi_{0}=0$. Therefore, the phase difference $\psi$ can be interpreted as the position of an overdamped particle moving in the tilted washboard potential $v(\psi)$ Lindner et al. (2004); Shlomovitz et al. (2014). The dynamic drift velocity $-\mathrm{d}v/\mathrm{d}\psi$ is symmetric around phase differences $\psi_{n}=(2n+1)\pi/2$ with $n\in\mathds{Z}$, which correspond to an order parameter of $\rho=1/2$. Therefore, if averaged over time, contributions from order parameters larger and smaller than $1/2$ exactly cancel out. In the case of finite phase discretizations $m$, the system is stochastic and it tends to spend a larger time in states with $\rho>1/2$. The reason for this can be understood by considering the Adler equation in the presence of noise and interpreting it as the governing equation of an overdamped Brownian particle in the potential $v(\psi)$. (In the case of the phase-discretized system, we may think of a ‘discrete’ potential whose increments determine the transition rates between states with different discrete phase differences.) For subcritical coupling strengths $\kappa<\kappa_{\mathrm{c}}$, the potential $v$ is (i) monotonic in $\psi$, (ii) convex in regions with $\rho>1/2$, and (iii) concave in regions with $\rho<1/2$; the latter can be seen by rewriting its second derivative as a function of the order parameter, $\mathrm{d}^{2}v/\mathrm{d}\psi^{2}=4\kappa(\rho-1/2)$. Therefore, the particle leaves regions with $\rho<1/2$ on the steepest slope of the potential, making it unlikely to return into the regions due to fluctuations whereas it leaves regions with $\rho>1/2$ where the potential is most shallow, rendering return events due to fluctuations much more likely.

IV.3 Oscillator precision at steady state

It is well-known that besides promoting synchronization, coupling can lead to an improvement of the oscillator precision, i.e., often damps frequency fluctuations Cross (2012). However, in the phase-discretized system, coupling not only tends to synchronize oscillators but is itself also a source of noise (cf. Eq. (15) and the discussion below). Hence, the effects of coupling on oscillator precision are not immediately obvious. To quantitatively assess these effects, we again consider the quality factor of the oscillators (see Section II), now for the coupled case: from the numerically obtained autocorrelation functions $G_{i}(t)=\langle x_{i}(t)x_{i}^{*}(0)\rangle$ of the two oscillators $i=1,2$, we obtain the quality factors by a fit with the exponential $\exp(\mathrm{i}\tilde{\omega}_{i}t-k_{i}t)$ as $Q_{i}=\tilde{\omega}_{i}/2\pi k_{i}$. From this, we compute the mean quality factor $\mathscr{Q}=(Q_{1}+Q_{2})/2$ as a proxy for the quality of the coupled system. Fig. 6 shows the steady-state order parameter $R$ and the steady-state quality factor $\mathscr{Q}$ as a function of the phase discretization $m$ and the coupling strength $\kappa$ for the case of equal frequencies (Fig. 6a,b) and the case of unequal frequencies (Fig. 6c,d). Remarkably, while the order parameter $R$ follows the general trends studied in the previous section, the quality factor $\mathscr{Q}$ exhibits certain optima along the coupling strength axis. In the case of equal frequencies (Fig. 6a,b), for a given phase discretization, increasing the coupling strength beyond the optimal value contributes more noise to the system than coupling is reducing. The location of this optimum depends on the intrinsic frequencies of the oscillators and for detuned frequencies, we consequently find two optima along the coupling strength axis (Fig. 6d). It is also interesting to note that there is no obvious correlation between synchrony and precision along the coupling strength axis, so that a high degree of phase synchrony can indeed be accompanied by large frequency fluctuations.

V Synchronization of many oscillators

We now turn to the dynamics of systems with larger numbers of oscillators and choose the classical case of an all-to-all coupled system to illustrate their behavior. For a system without self-coupling, the corresponding normalized adjacency matrix is given by $c_{ij}=(N-1)^{-1}(1-\delta_{ij})$.

V.1 ‘Mean-field’ formulation of the all-to-all coupled system

For an all-to-all coupling topology, the original Kuramoto model with sinusoidal coupling function can be rewritten in such a way that each oscillator individually couples to the order parameter $\Psi$, also called the ‘mean-field’ Strogatz (2000). The same is possible for the phase-discretized stochastic system specified by Eqs. (8–11) and (16), which can be rewritten in the form111 The rewriting in the ‘mean-field’ form relies on the fact that for sinusoidal coupling functions, the coupling term factorizes into a term containing the phase of the reference oscillator and the sum over all neighboring oscillators, e.g., $\smash{\sum_{j}\sin(\phi_{i}-\phi_{j})}=\smash{\operatorname{Im}(\mathrm{e}^{\mathrm{i}\phi_{i}}\sum_{j}\mathrm{e}^{-\mathrm{i}\phi_{j}})}=N\operatorname{Im}(\mathrm{e}^{\mathrm{i}\phi_{i}}\Psi)$ where $\Psi$ is the order parameter Eq. (17). The same rewriting can be applied to the transition rates of the phase-discretized system specified by Eqs. (11) and (16), which enables the representation Eqs. (27) and (28).

[TABLE]

where $\Psi$ is the Kuramoto order parameter defined in Eq. (17) and the coupling operator $\smash{\hat{\Gamma}^{\mathrm{MF}}_{i}}$ is given by

[TABLE]

where we have introduced the notation $[x]_{\pm}\equiv(1\pm x)/2$. Note that this rewriting also drastically reduces the computational effort to simulate the model222For stochastic simulations, the advantage of the form given by Eqs. (27) and (28) is that in order to compute all reaction propensities, it is sufficient to compute the order parameter $\Psi$ and the quantity given by Eq. (28) for each of the $N$ phases instead of computing all $N(N-1)$ pairwise phase differences.

Likewise, the corresponding linear noise approximation Eqs. (12–15) can be recast in the form

[TABLE]

where $\mu_{i}=|\omega_{i}|+\kappa$ is the effective noise strength for oscillator $i$, $\eta_{i}$ is Gaussian white noise with $\langle\eta_{i}(t)\rangle=0$ and $\langle\eta_{i}(t)\eta_{j}(t^{\prime})\rangle=\delta_{ij}\delta(t-t^{\prime})$, and where we have used the definition of the two global quantities

[TABLE]

The first quantity is the Kuramoto order parameter associated with the ‘macroscopic’ phases $\Phi_{i}$ and the second one convolves the macroscopic phases with the random variables $\xi_{i}$.

V.2 Synchronization transient

As for the case of two coupled oscillators, we assess the synchronization transient and the steady-state phase coherence for the many-oscillator system. To this end, we consider the case of homogeneous frequencies $\omega_{i}=\omega$. Fig. 7a illustrates the synchronization transient by showing the time-dependent order parameter $\rho(t)=\langle|\Psi(t)|^{2}\rangle$ for different numbers of oscillators and phase discretizations (cf. Fig. 2a). Figs. 7b,c show the synchronization time $\tau_{\mathrm{h}}$, Eq. (20), to reach an order parameter of $\rho=1/2$ from complete desynchronization as well as the time $\tau_{\nu}$, Eq. (21), to reach a fraction $\nu=0.95$ of the steady-state order parameter $R$. Note that for coarse phase discretizations, the system may not reach an order parameter of $1/2$ at all, in which case the time $\tau_{\mathrm{h}}$ is undefined (Fig. 7b). Generally, the shown synchronization times, which characterize the nonlinear transient from complete desynchronization to synchrony, increase with the number of oscillators $N$ but can exhibit a nonmonotonic behavior in the phase discretization $m$ for a large enough number of oscillators.

V.3 Steady-state phase coherence and oscillator precision

Turning to the steady-state phase coherence, Fig. 7d shows the steady-state order parameter $R$ as a function of the phase discretization for different numbers of oscillators (dots) as well as a comparison with the linear noise approximation (curves), Eqs. (29–31). While the order parameter increases monotonically in the phase discretization, the behavior for larger numbers of oscillators hinting at a synchronization transition for a finite value of the phase discretization (dark dataset and arrowhead in Fig. 7d). This transition is likely related to the synchronization transition of the classical Kuramoto model in the presence of noise, where partial synchrony is enabled when the noise strength drops below a critical value that depends on the coupling strength Acebrón et al. (2005). However, while in our case, the phase discretization is clearly related to the effective noise strength, cf. Eq. (7), it also introduces other effects such as sampling of the coupling function at discrete readout points, which may alter the behavior of the system apart from introducing fluctuations.

In the spirit of Section IV.3, we next assess the quality factor as the dimensionless ratio of the oscillation time scale and the exponential decay rate of the autocorrelation function. Fig. 7e shows the steady-state quality factor $Q$ computed from the average autocorrelation from all oscillators. We find a massive increase in oscillator precision when the phase discretization reaches values for which also the onset of partial synchrony is observed, cf. Fig. 7d. As expected, the quality factor also increases with the number of oscillators, a behavior that is well-known for coupled phase oscillator systems in general Cross (2012). Figs. 7d,e also suggest that for both the steady-state order parameter as well as the quality factor, the LNA specified by Eqs. (29–31) provides an excellent approximation in the limit of fine phase discretizations.

V.4 Onset of synchronization for inhomogeneous frequencies

Finally, we illustrate the behavior of many phase-discrete oscillators with inhomogeneous frequencies in the all-to-all coupled system Eq. (27). To this end, we consider an ensemble of systems with quenched disorder, i.e., with intrinsic transition frequencies $\omega_{i}$ drawn from a specific distribution $f(\omega)$ but fixed for each realization of the system. This scenario is well-studied for the classical Kuramoto model with unimodal and symmetric distributions $f$: In the thermodynamic limit of infinitely many oscillators, $N\to\infty$, this system exhibits a second-order synchronization phase transition at the critical coupling strength $\kappa_{\mathrm{c}}=2/(\pi f(0))$ Kuramoto (1984b); Strogatz (2000). We here draw the transition frequencies $\omega_{i}$ from a Cauchy distribution centered around zero,

[TABLE]

so that for the classical Kuramoto model in the thermodynamic limit, the phase transition occurs at $\kappa_{\mathrm{c}}=2$. Fig. 8 shows the order parameter $R$ as a function of the coupling strength $\kappa$ for different phase discretizations $m$ and the Kuramoto limit $m\to\infty$. Phase discretization decreases the limiting amount of synchrony and in some cases, even completely prohibits partial synchronization where the classical model is able to partially synchronize (see $m=5$ curve in Fig. 8). Again, the LNA given by Eqs. (29–31) provides good agreement with the phase-discretized model in the limit of fine discretizations.

VI Discussion

In this paper, we have presented a stochastic generalization of the Kuramoto model with discretized phases and investigated its synchronization behavior as well as its frequency fluctuations. Remarkably, while the phase-discretized model converges towards the deterministic Kuramoto dynamics in the limit of a continuous phase, many key observables exhibit a non-monotonic behavior. This leads to optima in the steady-state synchrony and oscillator precision for finite phase discretizations, which can exceed the corresponding values of the deterministic Kuramoto model. These features arise from an interplay of different effects that are a consequence of the phase discretization such as discrete sampling of the coupling function and the inherent stochasticity of the coupling process.

The discretization schemes introduced here enable a straightforward implementation of coupled stochastic oscillations in a Markov chain setting and can be useful in coupling cyclic dynamics to mesoscopic systems. Such systems might include, e.g., chemical reaction networks Gillespie (1977) and stochastic models of cell fate dynamics Klein et al. (2008), where cyclic processes may effectively depict periodic extrinsic signals such as the cell cycle Weber et al. (2014), circadian rhythms Zwicker et al. (2010); Brown (2014), or periodic signaling activity Sugimoto et al. (2012). Moreover, it is straightforward to computationally generalize the phase-discretized model to the case of non-Markovian transitions between phase states that entail non-exponential waiting time distributions Boguñá et al. (2014).

Here we have only taken a glimpse at the phenomena that can arise when phase-discretization of Kuramoto oscillators is combined with stochastic dynamics. To illustrate their behavior, we have chosen the generic cases of two oscillators and many oscillators with all-to-all coupling; therefore, we could not address the spatiotemporal dynamics of spatially distributed systems such as those with short-range (e.g., nearest-neighbor) interactions, which may give rise to interesting patterning phenomena Escaff et al. (2014). Moreover, we have chosen a generic but contingent discretization scheme for the coupling process (see Appendix B). It will be interesting to unfold the dynamics of different model realizations and to apply the proposed discretization schemes to, e.g., Kuramoto oscillators with inertia Tanaka et al. (1997); Olmi et al. (2014); Jörg (2015) and excitable dynamics Lindner et al. (2004) as well as time-delayed coupling Schuster and Wagner (1989); Yeung and Strogatz (1999); Earl and Strogatz (2003); Jörg et al. (2014) and signal filtering Pollakis et al. (2014); Wetzel et al. (2017), which goes beyond the Markovian approach.

Acknowledgements.

I thank B. D. Simons for discussions and L. Wetzel, L. G. Morelli, and I. M. Lengyel for critical comments on the paper. I acknowledge the support of the Wellcome Trust (grant number 098357/Z/12/Z).

Appendix A Linear noise approximation of the phase-discretized model with coupling

In this Appendix, we derive the linear noise approximation Eqs. (12–15). To this end, we perform a system size expansion of the system specified by Eqs. (8–11) in the standard way van Kampen (2011). The phase discretization $m$ is a natural candidate for the system size $\Omega$ as large $m$ lead to a more continuous phase. For the oscillator system with phases $\boldsymbol{\varphi}=(\varphi_{1},\ldots,\varphi_{N})$, we define the ‘macroscopic’ phases ${\boldsymbol{\Phi}}$ (that follow deterministic dynamics) and the random components ${\boldsymbol{\xi}}$ through the relation ${\boldsymbol{\varphi}}=\Omega{\boldsymbol{\Phi}}+\smash{\sqrt{\Omega}}{\boldsymbol{\xi}}$ where $\Omega=\varepsilon^{-1}=m/2\pi$. Furthermore, we define the probability distribution ${W}$ for the random components as ${W}(\boldsymbol{\xi},t)=P(\boldsymbol{\varphi}(\boldsymbol{\xi}),t)$. The next steps consist in calculating the time evolution of ${W}$, expanding in powers of $\sqrt{\Omega}$ and comparing coefficients. The coefficients of $\sqrt{\Omega}$ yield the equation

[TABLE]

whereas the coefficients of $\Omega^{0}$ result in

[TABLE]

where $\Gamma^{\prime}$ is the derivative of the coupling function. Eq. (33) describes the deterministic evolution of the macroscopic phases $\Phi_{i}$, while Eq. (34) is a Fokker–Planck equation for the random components $\xi_{i}$. The correspondence between Fokker–Planck and Langevin stochastic differential equations van Kampen (2011) enables to immediately write Eqs. (12–15) from Eqs. (33) and (34). In the case of no coupling, $\kappa=0$, Eqs. (13) and (14) reduce to $\mathrm{d}\Phi_{i}/\mathrm{d}t=\omega_{i}$ and $\mathrm{d}\xi_{i}/\mathrm{d}t=\smash{\sqrt{|\omega_{i}|}}\eta_{i}(t)$, so that the linear noise approximation Eq. (7) derived for Eq. (2) follows from Eqs. (12–15) and $\omega_{i}\geq 0$.

Appendix B Alternative generalizations of coupling

In this Appendix, we schematically discuss different possibilities to generalize the coupling term in a phase-discretized setting. To this end, we consider the Kuramoto model Eq. (1) for the case of two coupled oscillators without self-coupling and $\omega_{1}=\omega_{2}=0$. Schematically, the time evolution of oscillator $i=1,2$ can now be written as $\mathrm{d}\phi_{i}/\mathrm{d}t={K}_{i}$, where ${K}_{i}$ represents the dynamical frequency contribution from its coupling term. For simplicity, here we neither address the dependence of the ${K}_{i}$ on the phases nor their implicit time-dependence; these do not add any qualitative features to our considerations. We now illustrate different possibilities to generalize the coupling term by considering different stochastic processes (denoted by A, B, and C) for two discrete random variables $\varphi_{1}$ and $\varphi_{2}$ which all have in common that their expectation values satisfy $\mathrm{d}\langle\varphi_{i}\rangle/\mathrm{d}t={K}_{i}$.

For case A, we introduce four non-negative rates $k_{1}^{+}$, $k_{1}^{-}$, $k_{2}^{+}$, and $k_{2}^{-}$, with the only constraint that they satisfy $k_{i}^{+}-k_{i}^{-}={K}_{i}$. The stochastic dynamics is defined by the master equation

[TABLE]

where $P_{\mathrm{A}}=P_{\mathrm{A}}(\varphi_{1},\varphi_{2},t)$ and the $\hat{\varphi}_{i}^{\pm}$ are ladder operators, as defined in Eq. (10). Eq. (35) describes a system in which the forward and backward processes $\varphi\to\varphi+1$ and $\varphi\to\varphi-1$ are independent for each oscillator, leading to four parallel processes with rates $k_{i}^{\pm}$. In this case, stochastic reactions do not conserve the total number $\varphi_{1}+\varphi_{2}$. The coupling type investigated in this paper, Eq. (11), follows this spirit.

The stochastic dynamics of case B is defined by

[TABLE]

Eq. (36) describes a process in which for each oscillator at each point in time, depending on the sign of ${K}_{i}$ either the process $\varphi\to\varphi+1$ or $\varphi\to\varphi-1$ can occur, as indicated by the Heaviside function $\Theta$. A coupling in this spirit only admits a positive or negative frequency contribution at each point in time and importantly has zero contribution to the stochastic dynamics if ${K}_{i}=0$. This is not the case for coupling type A, where $K_{i}=0$ only imposes $k_{i}^{+}=k_{i}^{-}$.

Case C is only possible if ${K}_{1}={K}=-{K}_{2}$; this is the case, e.g., for symmetric bidirectional coupling $c_{ij}=c_{ji}$ and an odd coupling function such as $\Gamma(\phi)=\sin\phi$. Here we introduce two non-negative rates $k^{+}$ and $k^{-}$ with the only constraint that they satisfy $k^{+}-k^{-}={K}$ and define the stochastic dynamics by

[TABLE]

Eq. (37) describes a process in which the forward process $\varphi\to\varphi+1$ in one oscillator is always accompanied by a backward process $\varphi\to\varphi-1$ in the other oscillator, leading to the ‘exchange of phase quanta’ between the two oscillators and strict conservation of the total number $\varphi_{1}+\varphi_{2}$. It is clear that such a coupling type only works for symmetric coupling as any coupling-induced reaction will affect both oscillators involved.

This list of generalizations is by no means exhaustive and only gives a flavor of the different types of implementations of the stochastic coupling process. For instance, additional possibilities arise from the differences in how oscillators might internally process the coupling signals from different oscillators, e.g., whether they are processed independently Wetzel et al. (2017) or first integrated and then processed as a whole Pollakis et al. (2014). The adequate formalization to describe a specific system depends on the physical implementation of the coupling process.

Appendix C Steady-state order parameter and cross correlation of the two-oscillator system

In this Appendix, we derive Eq. (24) for the steady-state expectation value of the order parameter for two coupled phase-discretized oscillators. First, we obtain a master equation for the discrete phase difference $\theta=\varphi_{1}-\varphi_{2}$ by using Eq. (8) for $N=2$ and $c_{ij}=1-\delta_{ij}$ and marginalizing over one of the discrete phase variables, $\tilde{P}(\theta,t)=\sum_{\varphi_{1}}P(\varphi_{1},\varphi_{1}-\theta,t)$. For simplicity, we only consider the case $\omega_{1}\geq 0$, $\omega_{2}\geq 0$; the other cases follow analogously. Hence, we obtain the master equation for $\tilde{P}$ as

[TABLE]

where $\hat{\theta}^{\pm}$ are the ladder operators for the phase difference and where we have used the same convention for $[\cdot]_{\pm}$ as in Eq. (28). Here, $\tilde{\Gamma}(\phi)=(\Gamma(\phi)-\Gamma(-\phi))/2$ is the odd part of the coupling function and for the Kuramoto–Sakaguchi-type coupling Eq. (16), this evaluates to $\smash{\tilde{\Gamma}}(\phi)=\cos(\phi_{0})\sin(\phi)$. Next, we define the steady-state expectation values ${X}_{n}=\langle\mathrm{e}^{\mathrm{i}n\varepsilon\theta}\rangle$ and using the master equation (38), we obtain their time evolution as

[TABLE]

The key observation is that from the definition of the ${X}_{n}$ and $\varepsilon=2\pi/m$, it follows that ${X}_{0}=1$ and ${X}_{m}=1$, so that the set of equations given by (39) constitutes a closed hierarchy for the ${X}_{n}$ with $0\leq n\leq m$. At steady state, $\mathrm{d}{X}_{n}/\mathrm{d}t=0$, this yields the following set of algebraic equations,

[TABLE]

where the $\lambda_{n}$ are defined in Eq. (26). Solving this hierarchy starting from $n=m-1$, each expectation value ${X}_{n}$ can be expressed in terms of the next lower expectation value ${X}_{n-1}$. It can be shown by induction that this leads to the generic form

[TABLE]

where the $\Lambda_{n}$ satisfy the nonlinear recurrence relation

[TABLE]

with initial condition $\Lambda_{m}=0$. The continued fraction Eq. (25) is the solution to this recurrence relation as is obvious from repeatedly inserting Eq. (42) into itself. Since the cross correlation is given by ${X}={X}_{1}$, its exact solution is obtained from Eq. (41) as

[TABLE]

and Eq. (24) for the order parameter follows from this via Eq. (18) as $R=(1+\operatorname{Re}X)/2$.

Since the solution given by Eq. (43) is somewhat opaque due to its combinatorial complexity, we here give explicit expressions for ${X}$ for small $m$,

[TABLE]

where

[TABLE]