A common observation is that large groups of oscillatory biological units often have the ability to synchronize. A paradigmatic model of such behavior is provided by the Kuramoto model, which achieves synchronization through coupling of the phase dynamics of individual oscillators, while each oscillator maintains a different constant inherent natural frequency. Here we consider the biologically likely possibility that the oscillatory units may be capable of enhancing their synchronization ability by adaptive frequency dynamics. We propose a simple augmentation of the Kuramoto model which does this. We also show that, by the use of a previously developed technique [Ott and Antonsen, Chaos 18, 037113 (2008)], it is possible to reduce the resulting dynamics to a lower dimensional system for the macroscopic evolution of the oscillator ensemble. By employing this reduction, we investigate the dynamics of our system, finding a characteristic hysteretic behavior and enhancement of the quality of the achieved synchronization.

Examples from nature where many biological units in a large group oscillate and achieve synchrony are rhythmic clapping of humans in an audience, croaking of frogs, flashing of certain species of fireflies, and quorum sensing of yeast and bacteria. In this paper, we present a simple model for synchronization in such cases via a combination of both frequency and phase synchronization dynamics. A key feature of our model is that it can be rigorously reduced to a low dimensional dynamical system, thus greatly facilitating analysis and understanding. Results include the findings of a characteristic hysteretical behavior and enhancement of the quality of synchronization. The latter suggests that, to the extent that synchronization serves a significant functional purpose for some organisms, evolution may, in some cases, tend to promote the development of frequency adaptation dynamics.

Motivated by biological considerations, Winfree1 pioneered the modeling of synchronization of many coupled oscillators by means of a phase oscillator description in which the state of each oscillator i is given solely by its oscillation phase θi. Subsequently, Kuramoto2 presented a particularly appealing simplified version of the Winfree model

d θ i ( t ) d t = ω i + k N j = 1 N sin [ θ j ( t ) θ i ( t ) ] ,
(1)

where ωi is the natural oscillation frequency of oscillator i, which is a fixed, time-independent parameter giving the frequency of oscillator i in the absence of coupling to other oscillators ( j = 1 , 2 , , N ) , k is the “coupling constant” specifying the strength of the influence of other oscillators j on oscillator i, and we are interested in the case in which N ≫ 1. See Refs. 3–5 for discussion and generalization of Eq. (1). Reflecting the natural situations where no two of the biological units are exactly alike, the natural frequencies ωi are random draws from a continuous unimodal distribution function with some characteristic width Δω. We note that the spread in ωi opposes phase synchronization; e.g., for k = 0, given an initial condition where the θi(0) are bunched together about a particular phase value, in the absence of coupling (k = 0), the spread in natural frequencies leads to relaxation to an unbunched state on a timescale of order Δ ω 1 . On the other hand, the oscillator interaction term [the term in Eq. (1) proportional to k] promotes synchronization, e.g., for k sufficiently large, starting from an initial condition with small bunching, the interaction term grows, and as it does so, it “pulls” the phase θi of oscillator i more and more strongly toward the center of the bunch, thus increasing the amount of bunching.

Examples where coupling between biological units leads to synchronization include humans in an audience seeking to clap in unison,6–8 synchronous flashing of certain species of fireflies,9,10 synchronous chirping of crickets,11 and quorum sensing in yeast,12 among others. To the extent that the synchronization in some of these examples can promote functional goals of these organisms (e.g., mating), it may be that evolution has led to the emergence of mechanisms, in addition to that of phase pulling [the second term in Eq. (1)], that enhance synchronization. Partly motivated by this consideration, we study one such mechanism, “frequency adaptation.” Thus, we consider a slight modification of Eq. (1)

d θ i ( t ) d t = ω ̂ i ( t ) + k N j = 1 N sin [ θ j ( t ) θ i ( t ) ] .
(2)

That is, we replace the time-independent term ωi in Eq. (1) by a dynamically evolving frequency variable ω ̂ i ( t ) , and we will propose a dynamical model for the evolution of ω ̂ i ( t ) ( i = 1 , 2 , , N ) such that these frequency variables tend to be drawn to a common value. This same frequency-adaptation strategy has been previously employed in Refs. 6 and 10. However, an important point of our paper will be that, due to our simpler (but still reasonable) choice for modeling the frequency adaptation dynamics of ω ̂ i ( t ) , it will be possible to reduce the dynamics of our model to that of a low dimensional system (via the technique of Refs. 13–15) which can be analyzed in a rather complete way. We also hope that our paper may facilitate and simplify the related future analyses of other situations.

In formulating a model for the evolution of ω ̂ i ( t ) in Eq. (2) we assume that oscillator i can adjust its dynamic frequency ω ̂ i ( t ) based on its “observations” of the states of all the other oscillators as reflected by their global macroscopic state given by the complex order parameter

R ( t ) = r ( t ) e i Ψ ( t ) = 1 N j = 1 N exp [ i θ j ( t ) ] ,
(3)

where r may be thought of as a characterization of the strength or “quality” of the macroscopic synchronization of the oscillator ensemble (r = 0 corresponds to complete lack of synchrony, 0 < r < 1 corresponds to partial synchrony, and r = 1 corresponds to perfect synchrony), while Ψ gives the macroscopic ensemble phase.

We take as our frequency-adaptation model

d ω ̂ i ( t ) d t = ν ( r ( t ) ) [ ω ̂ i ( t ) d Ψ ( t ) d t ] γ ( r ( t ) ) [ ω ̂ i ( t ) ω i ] .
(4)

The coefficient ν(r) is a relaxation rate characterizing the tendency for oscillator i to be drawn toward (i.e., to adapt to) the current global macroscopic frequency (dΨ/dt) of the oscillator ensemble, while the relaxation rate γ(r) characterizes the tendency for the dynamic frequency of oscillator i to be drawn toward the inherent natural frequency (denoted ωi) of oscillator i. Since the first tendency depends on the existence of, at least, some synchrony of the oscillator ensemble, we suppose that ν(r) satisfies

ν ( 0 ) = 0 , d ν ( r ) / d r > 0.

On the other hand, we assume γ(0) and γ(r) are positive. For example, later, we shall consider the example

ν ( r ) = ν r , γ = γ 0 ,
(5)

where ν and γ0 are independent of r. Note that in the absence of synchrony (r = 0), Eq. (4) has the stable steady state ω ̂ i ( t ) = ω i .

Thus, our model consists of Eqs. (2)–(4), where the frequency adaptation mechanism is given by the first term on the right-hand side of (4), and our model reduces to the previous Kuramoto model (1) for ν = 0, ω ̂ i ( 0 ) = ω i . In contrast to the Kuramoto model, the state of node i is now given by two dynamical variables [ ω ̂ i ( t ) and θi(t)] rather than one [θi(t)].

We now consider solutions of Eq. (4) for some general time dependence of R = reiΨ. Whatever this time dependence is, the solution of the initial value problem (4), with the initial condition ω ̂ i ( 0 ) , is of the form

ω ̂ i ( t ) = Ω ( t ) + ω i H ( t ) + η i J ( t ) ,
(6)

where η i = ω ̂ i ( 0 ) ω i and the three time functions, Ω(t), H(t), and J(t), are the solutions of the following initial value problems:

(7)
d Ω ( t ) d t + [ ν ( r ) + γ ( r ) ] Ω ( t ) = ν ( r ( t ) ) d Ψ ( t ) d t ,
(7a)
with

Ω ( 0 ) = 0 ,
(7b)

(8)
d H ( t ) d t + [ ν ( r ) + γ ( r ) ] H ( t ) = γ ( r ( t ) ) ,
(8a)
with

H ( 0 ) = 1 ,
(8b)

(9)
d J ( t ) d t + [ ν ( r ) + γ ( r ) ] J ( t ) = 0 ,
(9a)
with

J ( 0 ) = 1.
(9b)

Because ν ( r ) 0 , γ ( r ) > 0 , the solutions of Eqs. (8) and (9) satisfy H(t) > 0 and J(t) > 0 for all finite time, and J decays to zero as t. Since J decays to zero as t, the initial frequency deviations ω ̂ i ( 0 ) ω i , are eventually forgotten and have no effect on the system attractors. However, in a case with more than one attractor, the initial frequency deviations can determine the attractor that is approached by a given orbit (e.g., see the last paragraph of Sec. IV).

Letting N, we introduce the distribution function f(θ, ω, η, t), where

g ( ω , η ) = 1 2 π 0 2 π f ( θ , ω , η , t ) d θ ,

is the marginal distribution function of the oscillator natural frequencies ω and initial frequency deviations η = ω ̂ ω . In terms of the joint probability distribution g(ω, η), the distribution of intrinsic frequencies is G 1 ( ω ) = d η g ( ω , η ) , and the distribution of initial dynamic oscillator frequencies ω ̂ i ( 0 ) is G 2 ( ω ̂ ) = d ω d η g ( ω , η ) δ [ ω ̂ ( ω + η ) ] . From Eqs. (2) and (6) the evolution equation for f(θ, ω, η, t) is

f t + θ { [ Ω ( t ) + ω H ( t ) + η J ( t ) + k r sin ( Ψ θ ) ] f } = 0.
(10)

Also, expressing the order parameter in Eq. (3) for N in terms of f, we have

R ( t ) = r ( t ) e i Ψ ( t ) = + d ω + d η 0 2 π d θ e i θ f .
(11)

Equations (10), (11), and (7)–(9) now form a complete closed system for the evolution of our oscillator ensemble with frequency adaptation.

We now seek solutions for this system in the form specified by the ansatz of Ref. 13 

f ( θ , ω , η , t ) = g ( ω , η ) 2 π { 1 + [ n = 1 α n ( ω , η , t ) e i n θ + ( c . c . ) ] } ,
(12)

where (c.c.) denotes the complex conjugate of the summation term. Following Refs. 14 and 15 and noting that H(t) is always positive, we conclude that, if f does not initially satisfy Eq. (12), it subsequently converges to (12) (where by convergence we mean convergence in the sense of Refs. 13 and 14). Inserting Eq. (12) in Eq. (10) we obtain

α t + k 2 ( R α 2 R * ) + i ( Ω + ω H + η J ) α = 0.
(13)

We now specialize to Lorentzian distributions in ω and η

g ( ω , η ) = Δ ω Δ η π 2 1 ( ω ω 0 ) 2 + Δ ω 2 1 η 2 + Δ η 2 .
(14)

(This corresponds to a Lorentzian distribution for G 2 ( ω ̂ ) of width Δω + Δη.) We can perform the ω and η integrations in (11) by analytically continuing f into the complex ω and complex η planes, followed by closing of the ω and η integration paths with large semicircles along | ω | and | η | in the lower-half complex ω and η planes. This yields a contribution resulting from the poles of g at ω = ω0iΔω and η = –iΔη,

R * ( t ) = α ( ω 0 i Δ ω , i Δ η , t ) .
(15)

Setting ω = ω0iΔω and η = –iΔη in Eq. (13), it becomes

d R d t + k 2 ( | R | 2 1 ) R + [ i ( Ω + ω 0 H ) + ( Δ ω H + Δ η J ) ] R = 0 ,
(16)

or

d r / d t = ( k / 2 ) ( 1 r 2 ) r ( Δ ω H + Δ η J ) r ,
(17)
d Ψ / d t = Ω + ω 0 H .
(18)

Equations (17), (18) and (7)–(9) form a closed system for the macroscopic dynamics of the oscillator ensemble.

We now show that the solution of these equations leads to the phase Ψ of the ensemble increasing at the steady rate ω0, i.e., dΨ/dt = ω0 for all time. We first note from Eqs. (7b) and (8b), that Eq. (18) yields dΨ/dt = ω0 at t = 0. Further, from Eqs. (7a) and (8a), we have d ( Ω + ω 0 H ) / d t + ( ν + γ ) ( Ω + ω 0 H ) = ν ( d Ψ / d t ) + ω 0 γ . Then, using Eq. (18)

d 2 Ψ d t 2 + γ d Ψ d t = γ ω 0 .
(19)

The solution to Eq. (19) subject to the initial condition dΨ/dt = ω0 at t = 0 is dΨ/dt = ω0 for all time.

Setting d/dt = 0 in Eqs. (17), (7)–(9) and further specializing to the case where ν(r) and γ(r) are given by Eq. (5), Eq. (17) yields

[ k 2 ( 1 r 0 2 ) Δ ω γ 0 ν r 0 + γ 0 ] r 0 = 0 ,
(20)

where r0 denotes a steady fixed point for the dynamics of r(t). Figure 1 shows plots of the solutions of (20) shown as solid lines (stable solutions) and dashed lines (unstable solutions). Each of these plots is for a different value of ν / γ 0 , namely ν / γ 0 = 0 , 1 , 2 , 4 , 16 , as labeled in the figure. The incoherent state (r0 = 0) bifurcates at a value of k = kc = 2Δω independent of the value of ν . Letting k = k* denote the k value at which the dashed and solid curves merge, we see that for k* < k < kc frequency adaptation introduces hysteresis whereby, in this range, r0 can either be in a stable synchronized state (r0 > 0) or in the incoherent state. In addition, as compared with the original solution of the Kuramoto model without frequency adaptation ( ν / γ 0 = 0 in Fig. 1), we see that for k > kc the quality of the synchronization, as indicated by the size of r0, can be much enhanced.

FIG. 1.

Steady state solutions versus the coupling constant k with the solid and dashed lines indicating the stable and unstable solutions. The values of ν / γ 0 for each curve are indicated by the arrows.

FIG. 1.

Steady state solutions versus the coupling constant k with the solid and dashed lines indicating the stable and unstable solutions. The values of ν / γ 0 for each curve are indicated by the arrows.

Close modal

We now investigate the dynamics of the order parameter leading to a stable equilibrium state for k in the hysteretic range k* < k < kc. We assume for simplicity that the initial spread in frequency deviations is zero (Δη = 0). In this case, the dynamics reduce from the three dimensional state space (r, H, J) to the two dimensional state space (r, H). It follows from Eq. (17) that the rate of change of the order parameter will be positive or negative if the quantity H is smaller or larger than k ( 1 r 2 ) / ( 2 Δ ω ) . Similarly, from Eq. (8), the rate of change of H will be positive or negative when H is smaller or larger than γ 0 / ( ν r + γ 0 ) .

These two conditions divide the unit square (r, H) into five regions as shown in Fig. 2. The curve shown in blue is where dH/dt = 0, and the one shown in black is where dr/dt = 0. The two curves cross at the two nonzero equilibrium points. In each region, we place a red arrow indicating the direction of the flow based on the signs of dH/dt and dr/dt. Flow lines are directed horizontally when crossing the blue line and vertically when crossing the black line. Equilibria are shown in Fig. 2 as solid circles. From the arrows, we see that the equilibrium at large r, as well as the incoherent equilibrium at (r, H) = (0, 1), are attractors, while the equilibrium at intermediate r is a saddle point. The stable manifold of the saddle point divides the plane into regions that are either attracted to the larger equilibrium point or to r = 0.

FIG. 2.

The curves d H / d t = H ̇ = 0 (blue) and d r / d t = r ̇ = 0 (black). The (red) arrows indicate the direction of the flow. This figure is for k/kc = 0.75 and ν / γ 0 = 4 .

FIG. 2.

The curves d H / d t = H ̇ = 0 (blue) and d r / d t = r ̇ = 0 (black). The (red) arrows indicate the direction of the flow. This figure is for k/kc = 0.75 and ν / γ 0 = 4 .

Close modal

In general, the exact shape of the stable manifold of the saddle, as well as the other trajectories in the (r, H) plane, depend on the parameters ν / γ 0 , k / k c and γ0ω. One limit is that of slow adaptation ( γ 0 , ν Δ ω ) . In this case, the trajectory in the (r, H) plane consists of a first phase of evolution of the order parameter r at constant H to the black curve r ̇ = 0 . This would then be followed by slow evolution along this curve in the direction dictated by the sign of H ̇ . In this case, any initial point with H greater than the value corresponding to the unstable equilibrium point will eventually evolve to r = 0. Any point with H less than the value of H for the unstable equilibrium will evolve to the stable equilibrium with positive r0.

Examination of the directions of trajectories in the (r, H) plane for k* < k < kc (Fig. 2) indicates that initial conditions starting with small values of the order parameter r and H(0) = 1 [Eq. (8b)] will always return to the incoherent state (0, 1). More generally, however, this need not be so. To see this, we revisit the choice [Eq. (14)] for the distribution g(ω, η). For this distribution, the spread of the ω ̂ i ( t ) from (6) is (ΔωH + ΔηJ) which, for H = 1, is always larger than Δω. Thus, as the contribution ηiJ(t) in (6) decays in time, the incoherent state is always stable. This can be changed by allowing η to be anti-correlated with ω such that the spread in instantaneous frequencies is initially small, and the order parameter grows. It is then possible, if the evolutions of H and J are sufficiently slow compared with the growth of r, to reach the large amplitude subcritical equilibrium from initial states with small r and k < kc.

In this paper, we have proposed a simple model for frequency adaptation of populations of coupled oscillatory units. As compared with the Kuramoto and Winfree models, this leads to an additional dynamically evolving variable in the specification of the state of each of the individual units. The proposed model is meant to capture possible dynamics of real biological systems, while at the same time being amenable to analytical attack, particularly by use of the techniques of Refs. 13–15.

The main finding from the solution of our model is that the frequency adaptation mechanism enhances the synchronization ability of coupled oscillatory units and leads to hysteretic onset of synchronization. Hysteretic onset has also been found in the previous paper, Ref. 6, using a different frequency adaptation model. The finding of hysteresis in two different models of frequency adaptation suggests that, to some extent, this feature may be characteristic of the general frequency adaptation mechanism rather than of the particular manner in which it is modeled.

Referring to Fig. 1, we see that for cases with frequency adaptation ( ν / γ 0 > 0 in Fig. 1), r0 is larger than for the case without adaptation ( ν / γ 0 = 0 in Fig. 1), meaning that frequency adaptation promotes synchronization. Thus, it seems plausible that evolution might, in some cases, favor the development of frequency adaptation in organisms for which synchronization has a functional benefit.

Furthermore, if it is supposed that hysteresis, as predicted by our frequency adaptation model, applies to the nocturnal, synchronous flashing of fireflies or croaking of frogs, then as night fell, assuming an initial unsynchronized state, the turn-on of synchrony would occur in a sudden (or “hard”) manner. Similarly, as conditions changed (e.g., dawn approached), the turn-off of synchrony would also be sudden. This is in contrast to the consequences of an assumption of Kuramoto-model-type dynamics for which the turn-on and turn-off transitions would be predicted to be “soft.”

Another aspect of our work has to do with the facts that the Kuramoto and Winfree models have been extended and generalized in various directions for the treatment of a great many different situations, and that it has been demonstrated that the resulting models can be treated using the ansatz of Refs. 13–15 (e.g., for networks of excitable and firing neurons, see Refs. 16–19) Thus, we hope that our paper may indicate a particular type of modeling that may, in some cases, be capable of accounting for additional dynamical state variables of individual units in a coupled ensemble.

This work was supported by the U.S. Army Research Office through Grant No. W911NF-12–1-0101.

1.
2.
Y.
Kuramoto
,
International Symposium on Mathematical Problems in Theoretical Physics
, Lecture Notes Physics Vol.
39
(
Springer
,
Berlin
,
1975
), p.
420
.
4.
E.
Ott
,
Chaos in Dynamical Systems
, 2nd ed. (
Cambridge University Press
,
Cambridge
,
2002
), Chap. 6.
5.
J.
Acebron
,
L.
Bonilla
,
C.
Vincenti
,
F.
Ritort
, and
R.
Spigler
,
Rev. Mod. Phys.
77
,
137
(
2005
).
6.
D.
Taylor
,
E.
Ott
, and
J.
Restrepo
,
Phys. Rev. E
81
,
046214
(
2010
).
7.
Z.
Néda
,
E.
Ravasz
,
T.
Vicsek
,
Y.
Brechet
, and
A.
Barabasi
,
Phys. Rev. E
61
,
6987
(
2000
).
8.
D.
Xenides
,
D.
Vlachos
, and
T.
Simos
,
J. Stat. Mech.
2008
,
P07017
(
2008
).
9.
J.
Buck
and
E.
Buck
,
Science
159
,
1319
(
1968
).
10.
B.
Ermentrout
,
J. Math. Biol.
29
,
571
(
1991
).
12.
S.
DeMonte
,
E.
d'Ovido
,
S.
Dano
, and
P.
Sorenson
,
Proc. Natl. Acad. Sci. U.S.A.
104
,
18377
(
2007
).
13.
E.
Ott
and
T.
Antonsen
,
Chaos
18
,
037113
(
2008
).
14.
E.
Ott
and
T.
Antonsen
,
Chaos
19
,
023117
(
2009
).
15.
E.
Ott
,
B.
Hunt
, and
T.
Antonsen
,
Chaos
21
,
025112
(
2011
).
16.
T.
Luke
,
E.
Barreto
, and
P.
So
,
Neural Comput.
25
,
3207
(
2013
).
17.
D.
Pazó
and
E.
Montbrió
,
Phys. Rev. X
4
,
011009
(
2014
).
18.
C.
Laing
,
Phys. Rev. E
90
,
010901
(
2014
).
19.
S.
Chandra
,
D.
Hathcock
,
K.
Crain
,
T.
Antonsen
,
M.
Girvan
, and
E.
Ott
,
Chaos
27
,
033102
(
2017
).