Kuramoto Model with Delay: The Role of the Frequency Distribution

Abstract

The Kuramoto model is a classical model used for the describing of synchronization in populations of oscillatory units. In the present paper we study the Kuramoto model with delay with a focus on the distribution of the oscillators’ frequencies. We consider a series of rational distributions which allow us to reduce the population dynamics to a set of several delay differential equations. We use the bifurcation analysis of these equations to study the transition from the asynchronous to synchronous state. We demonstrate that the form of the frequency distribution may play a substantial role in synchronization. In particular, for Lorentzian distribution the delay prevents synchronization, while for other distributions the delay can facilitate synchronization.

1. Introduction

In many complex systems of various nature a similar pattern of collective behaviour can be observed: The adjustment of rhythms of oscillating systems due to their interaction. This phenomenon is called synchronization and is ubiquitous in the real world with examples ranging from the synchronous firing of pacemaker cells in the heart and neurons in the brain to the synchronous rotation of electric generators in power grids [1,2,3]. For a long time it was a mystery how the synchronization can emerge despite the inevitable diversity in the natural frequencies of different units. It was Kuramoto who introduced a mathematical model of coupled oscillators that allowed this problem to be solved [4,5,6]. Motivated by the behavior of chemical and biological oscillators, this model later turned out to be quite general and applicable to such systems as coupled arrays of Josephson junctions [7] or populations of of biological neurons [8,9].

What is now known as the “Kuramoto model” consists of N phase oscillators with the harmonic interaction function:

(1)$\dot{{\theta }_j}\left(t\right)={\omega }_j-\frac{K}{N}\sum _{k}sin\left({\theta }_j-{\theta }_k\right),$

where $j=1,\dots ,N$ is the unit number, ${\theta }_j\in S^1$ are the phase variables, ${\omega }_i$ are the natural frequencies, and K is the global coupling strength. In his pioneering work Kuramoto showed this model to be mathematically tractable and to demonstrate the phase transition from an asynchronous state to synchronization at the critical coupling strength.

In subsequent decades the Kuramoto model became a classical paradigm for studying synchronization phenomena and it found numerous applications (see reviews [10,11] for the examples). From the other side, the simplicity and generality of the Kuramoto model makes it perfect for various kinds of generalization and modification. One of the most natural ways to make the model more realistic is to include coupling delays which are inevitably in real world due to the finite speed of signal propagation [12,13,14]. To take into account coupling delays the phases ${\theta }_k\left(t\right)$ in the sum from the r.h.s. of (1) are replaced by their delayed versions ${\theta }_k(t-\tau )$. In the case of two coupled oscillators the delay causes multistabilty of synchronous solutions with different frequencies [15] (interestingly, multistability also emerges if the delay is introduced not into the coupling, but into the oscillators themselves [16]). The similar effect of multistability was found for globally coupled oscillators with time delays, which results in the possibility of discontinuous transitions between different regimes in addition to classical smooth transitions [17,18]. In the case of local coupling, the delays were shown to induce complex spatio-temporal patterns and clusters [19]. In the case of distant-dependant delays similar patterns emerge and become propagating structures [20,21]. The case of heterogeneous delays was also considered and the emergence of clusters with various phase relations was demonstrated [22]. In rings of oscillators with local coupling, the delay was shown to act differently depending on the network symmetry: It induces multistability for unidirectional coupling and stabilizes the most symmetric solution for bidirectional coupling [23]. In rings with nonlocal coupling, the delays give birth to travelling waves [24]. For networks of oscillators with complex connectivity, such as scale-free networks, coupling delays may induce explosive synchronization [25].

Another natural generalization of the Kuramoto model is the addition of external forces to the oscillator dynamics which results in the emergence of a phase-dependent and/or time-dependent term in the r.h.s. of (1). This term is often taken in the form $asin{\theta }_j$, in which case the phase oscillators turn into the so-called “active rotators”. The dynamics of active rotators is much richer than that of phase oscillators since the former can demonstrate either oscillatory or excitatory behaviour depending on the relation between the individual frequency ${\omega }_j$ and the parameter a which is often called a “non-isochronicity parameter”. Ensembles of globally coupled active rotators were first considered in Refs. [26,27], and two different collective regimes were found: Entrainment with the external force and mutual entrainment. The introduction of heterogeneous external forces leads to similar results [28]. The transition between locked and unlocked states can be similar to either a super-critical Andronov-Hopf bifurcation or a saddle-node bifurcation [29]. In the case of heterogeneous assemblies made up of excitable and oscillatory units, it has been demonstrated that the transition to synchrony may be classical or re-entrant, depending on the particular form of the frequency distribution [30]. In addition to the parameter regions with synchronous and asynchronous regimes, the regions of bistability between different regimes were found as well [31].

Other types of generalized Kuramoto models were studied as well, including time-dependent [8,32,33], adaptive [34,35,36,37,38] and memristive [39] coupling, including noise [26,27,40], non-harmonic coupling functions [41,42] and pulse coupling [43,44], systems on smooth manifolds [45,46,47], etc. Typically, the analytic study of the Kuramto model and its generalizations is performed in the thermodynamic limit when the number of units tends to infinity. In this case the microscopic equations (1) are replaced by the continuity equation, and the individual frequencies ${\omega }_j$ are replaced by the probability density $g\left(\omega \right)$. For the classical Kuramoto model the particular form of this distribution does not play any important role. First, the system is invariant under the change of variables ${\omega }_j\to {\omega }_j+\Omega$ for arbitrary $\Omega$, which means that the mean frequency $⟨\omega ⟩$ can be always set to zero. Second, for unimodal symmetric distributions, the critical coupling at which the transition takes place equals $K_c=2/\left(\pi g\left(0\right)\right)$, i.e., it depends only on the height of the distribution peak but not on its particular shape. However, when coupling delays or non-isochronicity are included, the shape of the distribution comes into play because it becomes important how the typical frequencies $\omega$ relate to the delay $\tau$ and non-isochronicity parameter a.

In the present paper we study how the synchronization transition in the Kuramoto model with delay depends on the shape of the frequency distribution $g\left(\omega \right)$. For this sake we consider a series of rational distributions which are all unimodal and symmetric but are different in the flatness of the peak and decay rate of the tails. Using the Ott-Antonsen approach [48,49] we derive low-dimensional dynamical systems governing the collective dynamics of the population and perform its bifurcation analysis. We show that the shape of the frequency distribution can play a significant role for the system properties. For widely used Lorentzian distribution, the delay always prevents synchronization by raising the critical coupling strength. However, for distributions with lighter tails, the delay can also promote synchronization by lowering the critical coupling.

2. Model

We consider a heterogeneous assembly of N oscillators with delayed coupling

(2)$\begin{array}{cc}\hfill \dot{{\theta }_j}\left(t\right)& ={\omega }_j-asin{\theta }_j\left(t\right)-\frac{K}{N}\sum _{k}sin\left({\theta }_j\left(t\right)-{\theta }_k(t-\tau )\right),\hfill \end{array}$

where $j=1,\dots ,N$ is the unit number, ${\theta }_j\in S^1$ are the phase variables, ${\omega }_j$ are the natural frequencies, a is the non-isochronicity parameter, K is the global coupling strength and $\tau$ is the coupling delay. Strictly speaking, ${\omega }_i$ can be called “natural frequencies” only for $a=0$. For non-zero a, the rotation becomes non-uniform with the frequency $\sqrt{{\omega }_j^2-a^2}$ for $0<a<|{\omega }_j|$. At $|{\omega }_j|=a$ an isolated unit undergoes a SNIPER bifurcation toward the excitatory regime. This very dependence of the local oscillatory dynamics on the parameter a allows us to call it a “non-isochronicity parameter”.

In order to characterize the degree of synchrony in the population we introduce the Kuramoto complex order parameter

(3)$R\left(t\right)=\frac{1}{N}\sum _j{e}^{i{\theta }_j\left(t\right)}.$

The absolute value of this parameter serves the main indicator of the system synchronization. When it is close to zero, the phases of the oscillators are not correlated, and the system is asynchronous. When the order parameter is sufficiently different from zero it indicates the emergence of a bunch of oscillators rotating with close phases, which means synchronization. Using the Kuramoto order parameter allows us to rewrite (2) as

(4)$\dot{{\theta }_j}={\omega }_j-\frac{a}{2i}(e^{i{\theta }_j}-e^{-i{\theta }_j})+\frac{K}{2i}(R_{\tau }{e}^{-i{\theta }_j}-R_{\tau }^{*}{e}^{i{\theta }_j}),$

where $R_{\tau }\equiv R(t-\tau )$ and the asterisk denotes the complex conjugate. In the thermodynamic limit $N\to \infty$, the macroscopic state of the system is described by the probability density function $f(\theta ,\omega ,t)$, which obeys the continuity equation

(5)$\frac{\partial f}{\partial t}+\frac{\partial }{\partial \theta }\left(fv\right)=0,$

with the velocity v being the r.h.s. of Equation (4). The Kuramoto parameter in this case is evaluated as

(6)$R\left(t\right)={\int }_{-\infty }^{\infty }d\omega {\int }_0^{2\pi }f(\theta ,\omega ,t)e^{i\theta }d\theta ,$

3. Reduction of the Collective Dynamics

Following the theory of Ott and Antonsen [48,49] we will look for the long-term dynamics of the continuity Equation (5) in the form

(7)$f(\theta ,\omega ,t)=\frac{g\left(\omega \right)}{2\pi }\left(1+\sum _{n=1}^{\infty }\left[{\left(z^{*}(\omega ,t)\right)}^n{e}^{in\theta }+{\left(z(\omega ,t)\right)}^n{e}^{-in\theta }\right]\right),$

where

(8)$g\left(\omega \right)={\int }_0^{2\pi }f(\theta ,\omega ,t)d\theta$

is the probability density function of the natural frequencies and

(9)$z(\omega ,t)={\int }_0^{2\pi }f(\theta ,\omega ,t)e^{i\theta }d\theta ,$

is the local complex order parameter of the subpopulation with the natural frequency $\omega$. Obviously, the global and the local order parameter are connected by the self-consistency condition

(10)$R={\int }_{-\infty }^{\infty }g\left(\omega \right)z\left(\omega \right)d\omega .$

Substituting (7) into (5), one obtains the following equations for $z(\omega ,t)$:

(11)$\dot{z}(\omega ,t)=i\omega z+\frac{a}{2}(1-z^2)+\frac{K}{2}(R_{\tau }-R_{\tau }^{*}{z}^2),$

which together with (10) defines a delay integro-differential equation describing the collective dynamics of the population in the thermodynamic limit.

As the next step we consider a family of rational distributions $g\left(\omega \right)$, namely

(12)$g_n\left(\omega \right)=\frac{c_n}{{(\omega -\Omega )}^{2n}+{\Delta }^{2n}},$

where n is natural, $\Omega$ is the mean frequency, $\Delta$ is the distribution half-width, and

(13)$c_n=\frac{1}{\pi }nsin\frac{\pi }{2n}{\Delta }^{2n-1}$

is the normalization constant. For $n=1$ this distribution turns into a classical Cauchy distribution, and for $n\to \infty$ it converges to a uniform distribution on the interval $\omega \in [\Omega -\Delta ;\Omega +\Delta ]$. Assuming the rational function $g\left(\omega \right)$ allows us to evaluate the integral (10) using the residue theorem. A similar approach was recently used for populations of quadratic integrate-and-fire neurons [50,51]. Consider the analytic extension of function $z(\omega ,t)$ to complex $\omega$, then the integration contour can be closed by an infinitely large arc in the upper complex half-plane. In this half-plane the function (12) has n simple poles

(14)${\omega }_k=\Omega +\Delta e^{i{\alpha }_k},$

where $k=\overline{1,n}$ and ${\alpha }_k=\pi (k-0.5)/n$. Thus, the integral (10) can be evaluated as

(15)$R\left(t\right)=-\frac{i}{\Delta }sin\frac{\pi }{2n}\sum _{k=1}^n({\omega }_k-\Omega )z({\omega }_k,t).$

Writing Equation (11) for ${\omega }_1$, ${\omega }_2$,…, ${\omega }_n$ allows us to obtain a closed set of n delay differential equations for complex variables

(16a)$\begin{array}{cc}\hfill \dot{z_k}& =i(\Omega +\Delta e^{i{\alpha }_k})z_k+\frac{a}{2}(1-z_k^2)+\frac{K}{2}(R_{\tau }-R_{\tau }^{*}{z}_k^2),\hfill \end{array}$

(16b)$\begin{array}{cc}\hfill R& =-isin\frac{\pi }{2n}\sum _{k=1}^n{e}^{i{\alpha }_k}{z}_k,\hfill \end{array}$

where $k=\overline{1,n}$ and $z_k\left(t\right)\equiv z({\omega }_k,t)$. Introducing $z_k=x_k+i{y}_k$ and $R=X+iY$ allows us to rewrite these equations in the real form:

$\begin{array}{cc}\hfill \dot{x_k}& =-\Omega y_k-\Delta (y_{k}cos{\alpha }_k+x_{k}sin{\alpha }_k)+\dots \hfill \end{array}$

(17a)$\begin{array}{cc}\hfill & \frac{a}{2}(1+y_k^2-x_k^2)+\frac{K}{2}(X_{\tau }(1+y_k^2-x_k^2)-2Y_{\tau }{x}_k{y}_k),\hfill \\ \hfill \dot{y_k}& =\Omega x_k+\Delta (x_{k}cos{\alpha }_k-y_{k}sin{\alpha }_k)-a{x}_k{y}_k+\dots \hfill \end{array}$

(17b)$\begin{array}{cc}\hfill & \frac{K}{2}(Y_{\tau }(1+x_k^2-y_k^2)-2X_{\tau }{x}_k{y}_k),\hfill \end{array}$

(17c)$\begin{array}{cc}\hfill X& =sin\frac{\pi }{2n}\sum _{k=1}^n(x_{k}sin{\alpha }_k+y_{k}cos{\alpha }_k),\hfill \end{array}$

(17d)$\begin{array}{cc}\hfill Y& =sin\frac{\pi }{2n}\sum _{k=1}^n(y_{k}sin{\alpha }_k-x_{k}cos{\alpha }_k).\hfill \end{array}$

This set of DDEs governs the collective dynamics of the population in the thermodynamic limit $N\to \infty$. The following analysis is based on system (17).

4. Studying the Role of the Coupling Delay

For the case of isochronous oscillations with $a=0$, the system always has a trivial steady state $z_k=R=0$ corresponding to asynchronous dynamics of the oscillators. This state is stable for weak coupling K and destabilizes via an Andronov-Hopf bifurcation when the coupling becomes sufficiently strong, which constitutes a classical Kuramoto scenario [5,10]. The stable limit cycle born in this bifurcation corresponds to a partial synchronization of the oscillators.

The dynamics of the system is illustrated in Figure 1 where the dynamics of the individual phases and the Kuramoto order parameter are shown for the synchronous and asynchronous regimes. In both cases, the system starts from random initial conditions. In the case of asynchronous dynamics, all the phases rotate incoherently, and the order parameter remains close to zero. When the synchronization is achieved, a bunch of oscillators quickly emerge whose phases rotate with the same frequency, and the order parameter reaches a sufficiently non-zero value.

In order to illustrate the transition from the asynchronous to the synchronous state we plot the dependence of the Kuramoto order parameter on the coupling strength in Figure 2. The order parameter is small for weak coupling and rapidly grows as the coupling strength exceeds the critical value. Note that the results obtained by the simulation of the reduced model (17) coincide with those obtained for the microscopic system up to the fluctuations induced by the finite size effects. Note also that adding of the coupling delay might sufficiently influence the system dynamics and shift the critical value of the coupling strength. Further, we will analyze the role of the delay in detail with the help of the reduced system.

In the case of zero delay $\tau =0$ the critical coupling $K_c$ depends on the distribution width $\Delta$ in a linear way. Indeed, according to the results of Kuramoto [52],

(18)$K_c=\frac{2}{\pi g(\Omega )}=\frac{2\Delta }{nsin\frac{\pi }{2n}}.$

For nonzero coupling delays the critical coupling becomes delay-dependent. In order to determine the bifurcation point it is necessary to write the characteristic equation for the trivial steady state which has the form $\left|D\right(\lambda \left)\right|=0$, where D is $\left(2n\right)\times \left(2n\right)$ matrix

(19)$D=\left(\begin{array}{cc}{D}^{xx}& D^{xy}\\ D^{yx}& D^{yy}\end{array}\right),$

and $D^{xx}$, $D^{xy}$, $D^{yx}$, $D^{yy}$ represent $n\times n$ matrices with the following elements:

(20)$\begin{array}{ccc}\hfill D_{km}^{xx}& =& (-\lambda -\Delta sin{\alpha }_k){\delta }_{km}+\frac{K}{2}sin\frac{\pi }{2n}{e}^{-\lambda \tau }sin{\alpha }_k,\hfill \end{array}$

(21)$\begin{array}{ccc}\hfill D_{km}^{xy}& =& (-\Omega -\Delta cos{\alpha }_k){\delta }_{km}+\frac{K}{2}sin\frac{\pi }{2n}{e}^{-\lambda \tau }cos{\alpha }_k,\hfill \end{array}$

(22)$\begin{array}{ccc}\hfill D_{km}^{yx}& =& (\Omega +\Delta cos{\alpha }_k){\delta }_{km}-\frac{K}{2}sin\frac{\pi }{2n}{e}^{-\lambda \tau }cos{\alpha }_k,\hfill \end{array}$

(23)$\begin{array}{ccc}\hfill D_{km}^{yy}& =& (-\lambda -\Delta sin{\alpha }_k){\delta }_{km}+\frac{K}{2}sin\frac{\pi }{2n}{e}^{-\lambda \tau }sin{\alpha }_k,\hfill \end{array}$

where ${\delta }_{km}$ equals one for $k=m$ and zero in other case.

At the Andronov-Hopf point a pair of roots $\lambda =\pm i\omega$ emerge, which allows us to determine the value of the coupling strength $K_b$ at the bifurcation point by solving $\left|D\right(i\omega \left)\right|=0$. For small delays, this equation can be solved numerically by taking (18) as the initial point, then the solution can be traced along the delay value as a parameter. The obtained dependence $K_b\left(\tau \right)$ is plotted in Figure 3a for the Lorentzian distribution of the oscillator frequencies ($n=1$). The bifurcation coupling shows a minimum at $\tau =0$ and grows rapidly and monotonically for non-zero delays. Note that we have calculated the bifurcation curve for both positive and negative delays. Although negative time delays are not physical, we use them in the bifurcation analysis for a reason that will become clear later. Namely, they will help us to find other bifurcation curves existing for positive delays.

An important point is that the existence of one Andronov-Hopf bifurcation curve in a system with time delay implies the existence of other bifurcation curves at other delay values due to the so-called reappearance of periodic solutions [53]. Indeed, if the bifurcation takes place at the coupling strength $K_b$ for the delay ${\tau }_0$ this means the existence of a limit cycle with the period $T=2\pi /\omega$ and vanishing amplitude. This implies the existence of the same limit cycle at the same coupling strength for the delays ${\tau }_k={\tau }_0+kT$, where k is an integer, which means that each of these points are also Andronov-Hopf points. Note, however, that the stability of the emergent limit cycle can change, which means that the bifurcation can be either supercritical or subcritical.

The bifurcation curve found by starting from the delay-less case reappears in the manner described above and leads to the emergence of other bifurcation curves as shown in Figure 3b. These curves demonstrate minimums at $\tau =k{T}_0$, where $T_0=2\pi /{\omega }_0$ and ${\omega }_0$ is the frequency at which the collective oscillations emerge for $\tau =0$. Obviously, it is the frequency of the distribution peak, i.e., ${\omega }_0=\Omega$. For non-zero delays the frequency becomes different, therefore the different bifurcation curves do not exactly match each other in shape. However, they all still a single minimum at the multiples of $T_0$. The trivial steady state, i.e., the asynchronous regime is stable when the coupling strength is below all the bifurcation curves. Thus, the synchronization border is defined by the lowest point of the curve and has a saw-like shape. This border coincides completely with that obtained in Ref. [17] for the same setting which corroborates the validity of our analysis.

The obtained results suggest that introduction of the coupling delay prevents the system synchronization: For non-zero delays, the critical coupling at which the oscillators start to synchronize increases with respect to the delay-free case (at best, the critical coupling does not change if the delay is a multiple of $T_0$). This result seems to be obvious from the physical viewpoint: It is harder to adjust if one receives outdated information. However, it turns out that the coupling delay can in some cases promote synchronization. This surprising effect is observed when the distribution $g\left(\eta \right)$ is different from Lorentzian, i.e., $n>1$.

In order to illustrate this effect we calculated the synchronization borders on the $\tau -K$ plane for different values of n. We adjust the distribution half-width as

(24)$\Delta =\frac{n}{2}sin\frac{\pi }{2n},$

so that the critical coupling for the delay-free case equals unity for all n. The results are plotted in Figure 4 and show a significant difference between the Lorentzian distribution and distributions with $n>1$. For the Lorentzian distribution the delay always prevents synchronization and the critical coupling is always not less than unity. For the distributions of higher order $n>1$, some delay values can promote synchronization so that the critical coupling becomes less than unity. Another feature of high-order distributions is the complex form of the synchronization border with many peaks and valleys of different shapes.

5. Conclusions

In this paper we performed an analysis of the Kuramoto model with coupling delay paying special attention to the distribution of the oscillator frequencies $\omega$. We used the method of Ott and Antonsen which allows one to reduce the system dynamics in the case of infinitely many oscillators. For the rational frequency distributions, the dynamics can be reduced to a set of delay-differential equations whose number equals the degree of the denominator. By the means of bifurcation analysis, we obtained the Andronov-Hopf bifurcation curves indicating a synchronization transition in the population and so constructed the synchronization border in the parameter plane. Our results have revealed the different role of the delay for different frequency distributions. Thus, for the Lorentzian distribution, the delay always prevents synchronization by increasing the critical coupling strength. In contrast, for the distributions different from Lorentzian, the delay can promote synchronization: For certain delay values, the critical coupling turns out to be lower than in the delay-less case.

In the studies of collective dynamics of heterogeneous populations, it is typical to consider a Lorentzian distribution of local parameters. The reason for this choice is the simplicity of analytical treatment: For example, in our study the Lorentzian distribution with $n=1$ leads to the reduced system (9) of just a single differential equation for the complex variable $z_1$. The Lorentzian distribution is often treated as paradigmatic and qualitatively reflects the properties of an arbitrary unimodal distribution. However, our results show that the particular shape of the distribution can play a significant role in the system behaviour and synchronization. In particular, the role of the coupling delay turns out to be opposite for the Lorentzian and non-Lorentzian distributions.

In the end, we would like to emphasize that the present study is limited to the local stability analysis and does not consider the global stability of the asynchronous state. It means that the stable asynchronous state might coexist with a synchronous state in some parameter regions leading to bistability areas, as was demonstrated in Ref. [17] for the Lorentz frequency distribution. The emergence of bistability is associated with the subcritical Andronov-Hopf bifurcation, while the supercritical bifurcation supports monostability. The type of the bifurcation can be determined by the calculation of the first Lyapunov coefficient [54,55] which could be one of the directions of the further investigation. Other possibilities include consideration of a broader class of frequency distributions, including non-unimodal ones.

Author Contributions

Conceptualization, V.V.K.; Methodology, V.V.K.; Software, A.A.Z.; Formal analysis, V.V.K. and A.A.Z.; Investigation, V.V.K. and A.A.Z.; Writing—original draft, V.V.K.; Supervision, V.V.K.; Funding acquisition, V.V.K. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The data that support the findings of this study are available on request from the corresponding author.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures

Figure 1. The dynamics of the system in the asynchronous (a,b) and synchronous (c,d) regimes. The top panels show the time traces of 10 randomly chosen phases [Forumla omitted. See PDF.], while the bottom panels show the time trace of the Kuramoto order parameter [Forumla omitted. See PDF.]. The coupling strength [Forumla omitted. See PDF.] for (a,b) and [Forumla omitted. See PDF.] for (c,d) The other parameters are [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.].

Figure 2. The dependence of the Kuramoto order parameter on the coupling strength. Black circles: [Forumla omitted. See PDF.], gray squares: [Forumla omitted. See PDF.]. The other parameters: [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.]. Solid lines indicate the results obtained by the simulation of the reduced system (17).

Figure 3. (a) The Andronov-Hopf bifurcation curve for system (17) with [Forumla omitted. See PDF.], [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] and [Forumla omitted. See PDF.]. (b) The same curve (black solid line) and its reappearing instances (black dashed lines). The gray thick line shows the synchronization border.

Figure 4. (a) The synchronization borders of system (17) for [Forumla omitted. See PDF.] (dotted line), [Forumla omitted. See PDF.] (dash-dotted line) and [Forumla omitted. See PDF.] (solid line). The mean frequency [Forumla omitted. See PDF.], the half-width [Forumla omitted. See PDF.] (b) Enlarged part of the panel (a).

References

1. Pikovsky, A.; Kurths, J.; Rosenblum, M.; Kurths, J. Synchronization: A Universal Concept in Nonlinear Sciences; Cambridge University Press: Cambridge, UK, 2003; Volume 12.

2. Strogatz, S.H. Sync: How Order Emerges from Chaos in the Universe, Nature, and Daily Life; Hachette: London, UK, 2012.

3. Winfree, A.T. The Geometry of Biological Time; Springer Science & Business Media: Berlin, Germany, 2013; Volume 12.

4. Kuramoto, Y. Self-entrainment of a population of coupled non-linear oscillators. Proceedings of the International Symposium on Mathematical Problems in Theoretical Physics; Kyoto, Japan, 23–29 January 1975; Springer: Berlin/Heidelberg, Germany, 1975; pp. 420-422.

5. Kuramoto, Y. Chemical Oscillations, Waves, and Turbulence; Springer: Berlin, Germany, New York, NY, USA, 1984; Volume 19, 156.

6. Kuramoto, Y.; Nishikawa, I. Statistical macrodynamics of large dynamical systems. Case of a phase transition in oscillator communities. J. Stat. Phys.; 1987; 49, pp. 569-605. [DOI: https://dx.doi.org/10.1007/BF01009349]

7. Wiesenfeld, K.; Colet, P.; Strogatz, S.H. Frequency locking in Josephson arrays: Connection with the Kuramoto model. Phys. Rev. E; 1998; 57, 1563. [DOI: https://dx.doi.org/10.1103/PhysRevE.57.1563]

8. Cumin, D.; Unsworth, C. Generalising the Kuramoto model for the study of neuronal synchronisation in the brain. Phys. D Nonlinear Phenom.; 2007; 226, pp. 181-196. [DOI: https://dx.doi.org/10.1016/j.physd.2006.12.004]

9. Breakspear, M.; Heitmann, S.; Daffertshofer, A. Generative models of cortical oscillations: Neurobiological implications of the Kuramoto model. Front. Hum. Neurosci.; 2010; 4, 190. [DOI: https://dx.doi.org/10.3389/fnhum.2010.00190] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/21151358]

10. Acebrón, J.A.; Bonilla, L.L.; Vicente, C.J.; Ritort, F.; Spigler, R. The Kuramoto model: A simple paradigm for synchronization phenomena. Rev. Mod. Phys.; 2005; 77, pp. 137-185. [DOI: https://dx.doi.org/10.1103/RevModPhys.77.137]

11. Rodrigues, F.A.; Peron, T.K.D.M.; Ji, P.; Kurths, J. The Kuramoto model in complex networks. Phys. Rep.; 2016; 610, pp. 1-98. [DOI: https://dx.doi.org/10.1016/j.physrep.2015.10.008]

12. Smith, H. An Introduction to Delay Differential Equations with Applications to the Life Sciences; Springer: Berlin/Heidelberg, Germany, 2011.

13. Erneux, T. Applied Delay Differential Equations; Surveys and Tutorials in the Applied Mathematical Sciences Springer Science & Business Media: Berlin/Heidelberg, Germany, 2009; Volume 3, 204.

14. Kashchenko, S.; Maiorov, V.V. Models of Wave Memory; Springer: Berlin/Heidelberg, Germany, 2015.

15. Schuster, H.G.; Wagner, P. Mutual Entrainment of Two Limit Cycle Oscillators with Time Delayed Coupling. Prog. Theor. Phys.; 1989; 81, pp. 939-945. [DOI: https://dx.doi.org/10.1143/PTP.81.939]

16. Kashchenko, A.A. Multistability in a system of two coupled oscillators with delayed feedback. J. Differ. Eqs.; 2019; 266, pp. 562-579. [DOI: https://dx.doi.org/10.1016/j.jde.2018.07.050]

17. Yeung, M.K.S.; Strogatz, S.H. Time Delay in the Kuramoto Model of Coupled Oscillators. Phys. Rev. Lett.; 1999; 82, pp. 648-651. [DOI: https://dx.doi.org/10.1103/PhysRevLett.82.648]

18. Choi, M.Y.; Kim, H.J.; Kim, D.; Hong, H. Synchronization in a system of globally coupled oscillators with time delay. Phys. Rev. E-Stat. Phys. Plasmas Fluids Relat. Interdiscip. Top.; 2000; 61, pp. 371-381. [DOI: https://dx.doi.org/10.1103/PhysRevE.61.371]

19. Nakamura, Y.; Tominaga, F.; Munakata, T. Clustering behavior of time-delayed nearest-neighbor coupled oscillators. Phys. Rev. E Stat. Phys. Plasmas Fluids Relat. Interdiscip. Top.; 1994; 49, pp. 4849-4856. [DOI: https://dx.doi.org/10.1103/PhysRevE.49.4849]

20. Zanette, D.H. Propagating structures in globally coupled systems with time delays. Phys. Rev. E-Stat. Phys. Plasmas Fluids Relat. Interdiscip. Top.; 2000; 62, pp. 3167-3172. [DOI: https://dx.doi.org/10.1103/PhysRevE.62.3167] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/11088811]

21. Jeong, S.O.; Ko, T.W.; Moon, H.T. Time-delayed spatial patterns in a two-dimensional array of coupled oscillators. Phys. Rev. Lett.; 2002; 89, 154104. [DOI: https://dx.doi.org/10.1103/PhysRevLett.89.154104] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/12365992]

22. Petkoski, S.; Spiegler, A.; Proix, T.; Aram, P.; Temprado, J.J.; Jirsa, V.K. Heterogeneity of time delays determines synchronization of coupled oscillators. Phys. Rev. E; 2016; 94, 12209. [DOI: https://dx.doi.org/10.1103/PhysRevE.94.012209] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27575125]

23. D’Huys, O.; Vicente, R.; Erneux, T.; Danckaert, J.; Fischer, I. Synchronization properties of network motifs: Influence of coupling delay and symmetry. Chaos; 2008; 18, 37116. [DOI: https://dx.doi.org/10.1063/1.2953582] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/19045490]

24. Laing, C.R. Travelling waves in arrays of delay-coupled phase oscillators. Chaos Interdiscip. J. Nonlinear Sci.; 2016; 26, 94802. [DOI: https://dx.doi.org/10.1063/1.4953663]

25. Peron, T.K.D.; Rodrigues, F.A. Explosive synchronization enhanced by time-delayed coupling. Phys. Rev. E; 2012; 86, 16102. [DOI: https://dx.doi.org/10.1103/PhysRevE.86.016102]

26. Shinomoto, S.; Kuramoto, Y. Phase Transitions in Active Rotator Systems. Prog. Theor. Phys.; 1986; 75, 1105. [DOI: https://dx.doi.org/10.1143/PTP.75.1105]

27. Sakaguchi, H.; Shinomoto, S.; Kuramoto, Y. Phase transitions and their bifurcation analysis in a large population of active rotators with mean-field coupling. Prog. Theor. Phys.; 1988; 79, pp. 600-607. [DOI: https://dx.doi.org/10.1143/PTP.79.600]

28. Arenas, A.; Vicente, C.J.P. Exact long-time behavior of a network of phase oscillators under random fields. Phys. Rev. E; 1994; 50, 949. [DOI: https://dx.doi.org/10.1103/PhysRevE.50.949]

29. Antonsen, T.M., Jr.; Faghih, R.T.; Girvan, M.; Ott, E.; Platig, J. External periodic driving of large systems of globally coupled phase oscillators. Chaos Interdiscip. J. Nonlinear Sci.; 2008; 18, 37112. [DOI: https://dx.doi.org/10.1063/1.2952447] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/19045486]

30. Lafuerza, L.F.; Colet, P.; Toral, R. Nonuniversal results induced by diversity distribution in coupled excitable systems. Phys. Rev. Lett.; 2010; 105, 084101. [DOI: https://dx.doi.org/10.1103/PhysRevLett.105.084101] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/20868099]

31. Klinshov, V.; Franovic, I.; Franović, I. Two scenarios for the onset and suppression of collective oscillations in heterogeneous populations of active rotators. Phys. Rev. E; 2019; 100, 062211. [DOI: https://dx.doi.org/10.1103/PhysRevE.100.062211] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/31962480]

32. Petkoski, S.; Stefanovska, A. Kuramoto model with time-varying parameters. Phys. Rev. E; 2012; 86, 46212. [DOI: https://dx.doi.org/10.1103/PhysRevE.86.046212]

33. Khatiwada, D.R. Numerical Solution of Finite Kuramoto Model with Time-Dependent Coupling Strength: Addressing Synchronization Events of Nature. Mathematics; 2022; 10, 3633. [DOI: https://dx.doi.org/10.3390/math10193633]

34. Timms, L.; English, L.Q. Synchronization in phase-coupled Kuramoto oscillator networks with axonal delay and synaptic plasticity. Phys. Rev. E; 2014; 89, 32906. [DOI: https://dx.doi.org/10.1103/PhysRevE.89.032906]

35. Ha, S.Y.; Noh, S.E.; Park, J. Synchronization of Kuramoto oscillators with adaptive couplings. SIAM J. Appl. Dyn. Syst.; 2016; 15, pp. 162-194. [DOI: https://dx.doi.org/10.1137/15M101484X]

36. Kasatkin, D.V.; Nekorkin, V.I. Dynamics of the Phase Oscillators with Plastic Couplings. Radiophys. Quantum Electron.; 2016; 58, pp. 877-891. [DOI: https://dx.doi.org/10.1007/s11141-016-9662-1]

37. Kasatkin, D.V.; Yanchuk, S.; Schöll, E.; Nekorkin, V.I. Self-organized emergence of multilayer structure and chimera states in dynamical networks with adaptive couplings. Phys. Rev. E; 2017; 96, 062211. [DOI: https://dx.doi.org/10.1103/PhysRevE.96.062211]

38. Feketa, P.; Schaum, A.; Meurer, T. Stability of cluster formations in adaptive Kuramoto networks. IFAC-Pap.; 2021; 54, pp. 14-19. [DOI: https://dx.doi.org/10.1016/j.ifacol.2021.06.141]

39. Xu, Q.; Liu, T.; Ding, S.; Bao, H.; Li, Z.; Chen, B. Extreme multistability and phase synchronization in a heterogeneous bi-neuron Rulkov network with memristive electromagnetic induction. Cogn. Neurodyn.; 2022; pp. 1-12. [DOI: https://dx.doi.org/10.1007/s11571-022-09866-3]

40. Zaks, M.A.; Neiman, A.B.; Feistel, S.; Schimansky-Geier, L. Noise-controlled oscillations and their bifurcations in coupled phase oscillators. Phys. Rev. E; 2003; 68, 66206. [DOI: https://dx.doi.org/10.1103/PhysRevE.68.066206]

41. Komarov, M.; Pikovsky, A. The Kuramoto model of coupled oscillators with a bi-harmonic coupling function. Phys. D Nonlinear Phenom.; 2014; 289, pp. 18-31. [DOI: https://dx.doi.org/10.1016/j.physd.2014.09.002]

42. O’Keeffe, K.P.; Strogatz, S.H. Dynamics of a population of oscillatory and excitable elements. Phys. Rev. E; 2016; 93, 62203. [DOI: https://dx.doi.org/10.1103/PhysRevE.93.062203] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27415251]

43. Canavier, C.C.; Tikidji-Hamburyan, R.A. Globally attracting synchrony in a network of oscillators with all-to-all inhibitory pulse coupling. Phys. Rev. E; 2017; 95, 032215. [DOI: https://dx.doi.org/10.1103/PhysRevE.95.032215]

44. Klinshov, V.; Lücken, L.; Yanchuk, S. Desynchronization by phase slip patterns in networks of pulse-coupled oscillators with delays: Desynchronization by phase slip patterns. Eur. Phys. J. Spec. Top.; 2018; 227, pp. 1117-1128. [DOI: https://dx.doi.org/10.1140/epjst/e2018-800073-7]

45. Fiori, S. A control-theoretic approach to the synchronization of second-order continuous-time dynamical systems on real connected Riemannian manifolds. SIAM J. Control Optim.; 2020; 58, pp. 787-813. [DOI: https://dx.doi.org/10.1137/18M1235727]

46. Markdahl, J.; Thunberg, J.; Goncalves, J. High-dimensional Kuramoto models on Stiefel manifolds synchronize complex networks almost globally. Automatica; 2020; 113, 108736. [DOI: https://dx.doi.org/10.1016/j.automatica.2019.108736]

47. Cafaro, A.D.; Fiori, S. Optimization of a control law to synchronize first-order dynamical systems on Riemannian manifolds by a transverse component. Discret. Contin. Dyn.-Syst.-Ser. B; 2022; 27, pp. 3947-3969. [DOI: https://dx.doi.org/10.3934/dcdsb.2021213]

48. Ott, E.; Antonsen, T.M. Low dimensional behavior of large systems of globally coupled oscillators. Chaos Interdiscip. J. Nonlinear Sci.; 2008; 18, 37113. [DOI: https://dx.doi.org/10.1063/1.2930766]

49. Ott, E.; Antonsen, T.M. Long time evolution of phase oscillator systems. Chaos Interdiscip. J. Nonlinear Sci.; 2009; 19, 23117. [DOI: https://dx.doi.org/10.1063/1.3136851] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/19566252]

50. Klinshov, V.; Kirillov, S.; Nekorkin, V. Reduction of the collective dynamics of neural populations with realistic forms of heterogeneity. Phys. Rev. E; 2021; 103, L040302. [DOI: https://dx.doi.org/10.1103/PhysRevE.103.L040302]

51. Pyragas, V.; Pyragas, K. Mean-field equations for neural populations with q -Gaussian heterogeneities. Phys. Rev. E; 2022; 105, 044402. [DOI: https://dx.doi.org/10.1103/PhysRevE.105.044402] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35590671]

52. Kuramoto, Y. Proceedings of the International symposium on Mathematical Problems in Theoretical Physics. Lect. Notes Phys.; 1975; 30, pp. 420-422.

53. Yanchuk, S.; Perlikowski, P. Delay and periodicity. Phys. Rev. E; 2009; 79, pp. 1-9. [DOI: https://dx.doi.org/10.1103/PhysRevE.79.046221] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/19518326]

54. Kuznetsov, Y.A. Elements of Applied Bifurcation Theory; Springer: New York, NY, USA, 2004.

55. Guckenheimer, J.; Holmes, P. Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields; Springer Science & Business Media: Berlin/Heidelberg, Germany, 2013; Volume 42.