Networks of coupled dynamical units give rise to collective dynamics such as the synchronization of oscillators or neurons in the brain. The ability of the network to adapt coupling strengths between units in accordance with their activity arises naturally in a variety of contexts, including neural plasticity in the brain, and adds an additional layer of complexity: the dynamics on the nodes influence the dynamics of the network and vice versa. We study a minimal model of Kuramoto phase oscillators including a general adaptive learning rule with three parameters (strength of adaptivity, adaptivity offset, adaptivity shift), mimicking learning paradigms based on spike-time-dependent plasticity. Importantly, the strength of adaptivity allows to tune the system away from the limit of the classical Kuramoto model, corresponding to stationary coupling strengths and no adaptation and, thus, to systematically study the impact of adaptivity on the collective dynamics. We carry out a detailed bifurcation analysis for the minimal model consisting of N = 2 oscillators. The non-adaptive Kuramoto model exhibits very simple dynamic behavior, drift, or frequency-locking; but once the strength of adaptivity exceeds a critical threshold non-trivial bifurcation structures unravel: A symmetric adaptation rule results in multi-stability and bifurcation scenarios, and an asymmetric adaptation rule generates even more intriguing and rich dynamics, including a period-doubling cascade to chaos as well as oscillations displaying features of both librations and rotations simultaneously. Generally, adaptation improves the synchronizability of the oscillators. Finally, we also numerically investigate a larger system consisting of N = 50 oscillators and compare the resulting dynamics with the case of N = 2 oscillators.

The synchronization of coupled oscillators is a fascinating manifestation of self-organization—indeed, self-emergent synchronization is a central process to a spectacular range of natural and technological systems, including the beating of the heart,16 flashing fireflies,17 pedestrians on a bridge locking their gait,18 genetic clocks,19 pendulum clocks hanging on a beam,20 mechanical oscillators,21,22 superconducting Josephson junctions,23 chemical oscillations,24,25 metabolic oscillations in yeast cells,26 life cycles of phytoplankton,27 and networks of neurons in the brain.3,4,28

A desirable property in real-world oscillator networks is the presence of synchronization whenever it ensures the proper functioning of the network: in the realm of technology, for instance, the AC between generators and consumers in a power grid needs to stay synchronized to ensure ideal power transmission and ultimately avoid power blackouts;5 and wireless networks require synchronized clocks29 to ensure safe data transmission; in the realm of biology, the proper functioning of the heart requires that the rhythmic electric activation of cardiac cells stays coordinated.16 

While the network interaction between dynamic units may give rise to intriguing collective behaviors such as synchronization, adding the ability to adapt the (coupling) weights on a network according to the dynamics on its nodes leads to co-evolutionary network dynamics.30 Indeed—“Intelligence is the ability to adapt to change”31—the adaptive dynamics of co-evolutionary networks may allow to increase their functional robustness.32,33 Adaptive co-evolutionary networks appear in a wide range of systems, including the vascular network,34,35 the glymphatic network of the brain,36,37 osteocyte network formation,38 social networks,39 and, in particular, in neural networks in the brain where neural plasticity plays an important role for learning,15,40 but also for the progression of certain neuro-degenerative diseases.41 In the current context, this raises the question if adaptivity in a coupled oscillator network can increase its ability to synchronize.

We devise a special variation of the Kuramoto model with adaptive coupling weights. While previous studies have considered such systems from various perspectives, we here introduce a parameterization that allows to systematically deviate from the traditional Kuramoto system with stationary uniform coupling and systematically study the effects of adaptation and the related bifurcation behavior.

This article is structured as follows. In Sec. II, we explain our adaptive Kuramoto model and how we parameterize it. In Sec. III, we consider N = 2 oscillators and carry out a detailed bifurcation analysis for several important parameter cases: (Sec. III A) Non-adaptive limit (classical Kuramoto model); (Sec.III B) Symmetric adaptation ( α = β = 0); (Sec. III C) Asymmetric adaptation. In Sec. IV, we numerically investigate larger oscillator systems with N = 50 and compare the resulting dynamic behavior with the results for N = 2. Section V concludes with a summary and discussion of our results.

A general model of N oscillators l [ N ] := { 1 , , N } with adaptive coupling is considered. The oscillator phase ϕ l = ϕ l ( t ) T = R / 2 π Z then evolves according to
d ϕ l d t = ω l + 1 N m = 1 N κ l m g ( ϕ m ϕ l ) ,
(1)
with intrinsic frequencies ω l. Oscillators l and m interact via the interaction function g and are coupled with the time-dependent coupling strength κ l m R. The coupling strengths evolve according to
d κ l m d t = ε ( A ( ϕ l ϕ m ) κ l m ) ,
(2)
where the adaptation or learning rule A = A ( ϕ ) defines how the coupling adapts according to the oscillator phases (states); i.e., the adaptation/learning rule is defined via a local interaction between oscillator’s phases. The second term with κ l m guarantees boundedness of the coupling strengths. The time scale at which adaptation may occur is set by ε: oscillator phases evolve on time scales 1 and coupling adapts on a time scale of 1 / ε.

1. Coupling function and adaptation rule

Specifically, we consider the simplest version of such a model and suppose that g and A are periodic functions in ϕ. We may think of these functions as truncations of Fourier series to first order, i.e., the coupling function becomes
g ( ϕ ) = sin ( ϕ + α ) ,
(3)
and the adaptation rule (or function) becomes
A ( ϕ ) = a 0 + a 1 cos ( ϕ + β ) ,
(4)
where α and β tune the level of cosine vs sine interaction in the coupling function and adaptation rule, respectively. Thus, the oscillator dynamics is defined by the Kuramoto–Sakaguchi model with phase-lag α; the adaptation/learning rule has an adaptation offset, a 0 0, and the adaptivity a 1 (i.e., the strength of adaptation). The adaptation rule is inspired by models of spike-time-dependent plasticity, a Hebbian synaptic learning rule,42 and can be approximated in phase oscillator networks43 as we do here. The shape of the adaptation rule induces different types of adaptation or learning.9,10,12 In this context, the adaptation shift β tunes the type of interaction. Restricting our attention to positive adaptivity, a 1 > 0, the following adaptation regimes may be distinguished (see also Fig. 1):
  1. Symmetric adaptation with in-phase gain ( β = 0) amplifies (or suppresses) the coupling strength between oscillators that have a phase difference close to 0 (or π).

  2. Symmetric adaptation with anti-phase gain ( β = π) amplifies (or suppresses) the coupling strength between oscillators that have a phase difference close to π (or 0).

  3. Non-symmetric adaptation ( β 0 , π) amplifies (or suppresses) coupling strengths between oscillators with non-trivial phase difference; e.g., for the odd adaptation function with β = π / 2, coupling strengths associated with phase difference close to π / 2 (or π / 2) are amplified (or suppressed).

FIG. 1.

The adaptation rules A ( ϕ ) for distinct values of adaptation shifts β (shown for a 1 > 0).

FIG. 1.

The adaptation rules A ( ϕ ) for distinct values of adaptation shifts β (shown for a 1 > 0).

Close modal

Note that the system has the parameter symmetry ( a , β ) ( a , β + π ). Thus, if we consider β = 0 for both positive ( a > 0) and negative ( a < 0) adaptivity, we need not include the case of β = π in our investigation. In addition to these special cases, we are also interested in investigating mixed-type learning rules such as the case β = π / 4, see Sec. III C 3.

Finally, note that symmetric/asymmetric adaptation corresponds to symmetric/asymmetric coupling, i.e., undirected/directed network weights.

2. Reduced/rescaled model and its symmetries

Furthermore, we note that the “self-couplings,” κ l l ( t ), are decoupled from the oscillator phases; vice versa, if g ( 0 ) = sin α = 0, oscillators are independent of the self-coupling. But even if sin α 0, the post-transient or asymptotic values κ l l ( ) a 0 + a 1 cos β can be absorbed into the intrinsic frequencies, ω l + sin α κ l l ( ) ω l. We may, therefore, disregard all coupling terms κ l l.44 Furthermore, we may reduce the number of parameters involved by appropriately rescaling parameters and variables,
a 0 t t , ω l a 0 ω l , κ l m a 0 κ l m , ε a 0 ε .
Since ε > 0, we restrict ourselves to a 0 > 0. Introducing the (relative) adaptivity, a := a 1 / a 0, we then obtain the governing equations for the model,
d ϕ l d t = ω l + 1 N m = 1 N κ l m sin ( ϕ m ϕ l + α ) ,
(5a)
d κ l m d t = ε ( 1 + a cos ( ϕ l ϕ m + β ) κ l m ) ,
(5b)
where l m.

Since g ( ϕ ) and A ( ϕ ) only depend on phase differences, (5) is invariant to shifts in constant phase and frequency, i.e., ϕ l ( t ) ϕ l ( t ) + Φ + Ω t with constants Φ , Ω R.

3. Non-adaptive case (a= 0)

When a = 0, the model asymptotically strives to κ l m = 1 for all l , m [ N ], corresponding to the Kuramoto model with (rescaled) uniform stationary coupling strength, κ = 1. The emergent synchronization in the model can be characterized by the order parameter
Z ( t ) = 1 N l = 1 N e i ϕ l ( t ) .
(6)
The Kuramoto model exhibits the following asymptotic behavior.7,32,45 In the continuum limit, N , when the coupling strength is subcritical ( κ < κ c), the order parameter converges to 0, | Z | 0, corresponding to incoherent oscillations; vice versa, for supercritical coupling ( κ > κ c), the order parameter approaches a non-zero value corresponding to partial synchronization, | Z | C > 0. (Note that for finite oscillator systems, the order parameter is subject to pseudo-random fluctuations and it is preferable to consider the time-averaged order parameter, | Z ( t ) | t. In the subcritical regime, fluctuations will prevent the order parameter from precisely attaining 0.) For the continuum limit N , the threshold for this transition occurs at κ c = 2 / ( π g ( 0 ) ),32,33 where g ( ω ) is a unimodal frequency distribution with maximum g ( 0 ) and width σ. Note that we normalized the coupling to 1; thus, the synchronization threshold, σ c, is indirectly determined via 1 = 2 / ( π g ( 0 ) ). Vice versa, if a 0, we allow the network to be adaptive and we move away from the manifold corresponding to the Kuramoto model with stationary coupling. We will also use the order parameter (6) to characterize the adaptive dynamics for N = 50 in Sec. IV.
We consider the minimal network with N = 2 oscillators. Introducing ϕ := ϕ 1 ϕ 2 and ω := ω 1 ω 2, the governing equations become
d ϕ d t = ω + 1 2 κ 12 sin ( α ϕ ) 1 2 κ 21 sin ( α + ϕ ) ,
(7a)
d κ 12 d t = ε ( 1 + a cos ( β + ϕ ) κ 12 ) ,
(7b)
d κ 21 d t = ε ( 1 + a cos ( β ϕ ) κ 21 ) .
(7c)
This system of equations has the additional symmetries
( a , β ) ( a , β + π ) ,
(8a)
( a , α , ϕ ) ( a , α + π , ϕ + π ) ,
(8b)
( α , β , κ 12 , κ 21 ) ( α , β , κ 21 , κ 12 ) .
(8c)
Sometimes, it is useful to consider the dynamics in terms of the difference Δ := κ 12 κ 21 and sum Σ := κ 12 + κ 21. Equations (7b) and (7c) then become
d Δ d t = ε ( a [ cos ( β + ϕ ) cos ( β ϕ ) ] Δ ) ,
(9a)
d Σ d t = ε ( 2 + a [ cos ( β + ϕ ) + cos ( β ϕ ) ] Σ ) .
(9b)

Due to (8a), we may restrict the parameter range to π / 2 < β π / 2; due to (8b), we may restrict the parameter range to π / 2 < α π / 2; furthermore, (8c) allows to restrict either the range of β or α to positive values, and we chose to restrict 0 β π / 2. However, observe that all symmetry transformations also affect other parameters and variables; e.g., the symmetry (8b) effectively swaps (near-)phase-locked states where ϕ is close to zero, with (near-)antiphase states where ϕ is close to π. Thus, if a phase-locked stationary state is stable for α = 0 with adaptivity a, then an antiphasic stationary state is also stable for α = π with a. There are several other symmetries, see the  Appendix, some of which are reflected in the stability diagrams [Figs. 3–6(a)].

FIG. 2.

Bifurcation diagrams (equilibria only) in ( ω , κ ) and ( ω , ϕ ) for varying values of a (indicated as green lines in Fig. 3) with α = β = 0 and ε = 0.2. The diagrams show stable nodes (black curves), saddles, unstable nodes (both black dashed curves), stable spirals (red lines), and unstable spirals (red dashed curves). Unstable limit cycles (librations; not shown) emerge from subcritical Hopf bifurcations. SN and HB denote saddle-node and Hopf points, respectively.

FIG. 2.

Bifurcation diagrams (equilibria only) in ( ω , κ ) and ( ω , ϕ ) for varying values of a (indicated as green lines in Fig. 3) with α = β = 0 and ε = 0.2. The diagrams show stable nodes (black curves), saddles, unstable nodes (both black dashed curves), stable spirals (red lines), and unstable spirals (red dashed curves). Unstable limit cycles (librations; not shown) emerge from subcritical Hopf bifurcations. SN and HB denote saddle-node and Hopf points, respectively.

Close modal
FIG. 3.

Stability diagram in ( ω , a ) for α = β = 0 with ε = 0.2. ( ω , a ) parameter regions are shaded according to which stable states they contain Drift (D) and/or locked states (L). BT and CP denote Bogdanov–Takens and Cusp codimension 2 bifurcation points. Het, heteroclinic bifurcation with saddle and drift cycle and D-loss stability loss of drift cycle.

FIG. 3.

Stability diagram in ( ω , a ) for α = β = 0 with ε = 0.2. ( ω , a ) parameter regions are shaded according to which stable states they contain Drift (D) and/or locked states (L). BT and CP denote Bogdanov–Takens and Cusp codimension 2 bifurcation points. Het, heteroclinic bifurcation with saddle and drift cycle and D-loss stability loss of drift cycle.

Close modal

Furthermore, we consider ε > 0. As our analysis will show, we find interesting nontrivial dynamic behavior (at least) if ε = 0.2. We use this value throughout the analysis unless specified otherwise.

Certain parameter choices lead to effectively lower-dimensional dynamics. It is instructive to first consider these cases (Secs. III A and III B) before analyzing the more general case leading to three-dimensional dynamics (Sec. III C).

Non-adaptive dynamics is obtained for two limiting cases: either ε = 0, so that κ 12 and κ 21 are constants; or the adaptivity is zero, a = 0, and the coupling strengths asymptotically become identical with κ 12 ( t ) , κ 21 ( t ) 1 as t , irrespective of the value of β. In either cases, the system is effectively one-dimensional with the dynamics given by (7a), i.e., the system corresponds to the classical Kuramoto model. Since we are only interested in ε > 0, we consider the case of a = 0 and let κ 12 = κ 21 = 1 to consider the post-transient dynamics. Equation (7a) can be cast as
d ϕ d t = ω cos α sin ϕ .
(10)
Since we restricted α ( π / 2 , π / 2 ], we have cos α 0. Excluding the effectively decoupled case with cos α 0, we may rescale time and ω with cos α to obtain
d ϕ d t = ω sin ϕ .
(11)
When the frequency mismatch exceeds the coupling strength, | ω | > 1, there are no fixed points and the distance ϕ = ϕ 1 ϕ 2 keeps increasing, amounting to drifting oscillators. When the frequency mismatch is smaller, | ω | < 1, we have the two equilibria
ϕ := arcsin ω , ϕ + := π arcsin ω ,
(12)
ϕ is a stable equilibrium while ϕ + is unstable (note that stability types are inverted for cos α < 0). These equilibria correspond to (frequency-)locked states with the two oscillators having constant phase differences. When the oscillators share their intrinsic frequencies, ω = 0, the oscillators are phase-locked. For | ω | = 1, the two equilibria collide in a saddle-node bifurcation on an invariant cycle (SNIC). These drift/locked states, as well as the saddle-node bifurcation, remain as a basic dynamic structure when the system is adaptive with a 0.

1. Reduction

The adaptation rule is symmetric (or an even function) for β = 0 and β = π, so that Eqs. (7b) and (7c) attain identical structure. As a result, the three-dimensional dynamics of (7) is asymptotically described by dynamics on a two-dimensional subspace with symmetric coupling, κ 12 = κ 21, corresponding to undirected network weights. To see this, we substitute β = 0 into Eq. (9), which simplifies to
d Δ d t = ε Δ ,
(13a)
d Σ d t = ε ( 2 [ 1 + a cos ϕ ] Σ ) .
(13b)
Defining κ ( t ) := κ 12 ( t ) = κ 21 ( t ) = Σ ( t ) / 2, (7a) and (13b) become
d ϕ d t = ω κ cos α sin ϕ ,
(14a)
d κ d t = ε ( 1 + a cos ϕ κ ) .
(14b)
This system describes the dynamics on a two-dimensional subspace, i.e., the plane defined by κ 12 = κ 21 in the full ( ϕ , κ 12 , κ 21 ) phase space. As is evident from inspecting (13a), this plane is globally attracting since Δ 0 as t .

Equation (14) admits for two degenerate cases: (i) if ω = cos α = 0, ϕ is constant and, thus, κ ( t ) 1 + a cos ϕ exponentially fast and (ii) if cos α = 0 with ω 0, we have ϕ ( t ) = ϕ 0 + ω t so that oscillators drift apart and κ ( t ) forever undergoes an oscillation.

Since we restricted α ( π / 2 , π / 2 ], we may assume cos α 0. Further restricting cos α 0, we may rescale time and related parameters with cos α to obtain
d ϕ d t = ω κ sin ϕ ,
(15a)
d κ d t = ε ( 1 + a cos ϕ κ ) .
(15b)

2. Stability analysis

Equilibrium conditions are given by Eqs. (15b) and (15a) with
κ = 1 + a cos ϕ
(16)
and
ω = κ sin ϕ = ( 1 + a cos ϕ ) sin ϕ .
(17)
Using Euler’s identity, we may eliminate ϕ to obtain an equilibrium condition in κ only,
1 = ω 2 κ 2 + ( κ 1 ) 2 a 2 .
(18)
The equilibria for ϕ are given as
tan ϕ = ω κ a κ 1 ,
(19)
as a function of κ. Equilibria (and their bifurcations) are shown in Fig. 2 for varying values of ω and a.
To determine the asymptotic stability of these equilibria, we consider the Jacobian of (15),
J = [ κ cos ϕ sin ϕ ε a sin ϕ ε ] .
(20)
Using equilibrium conditions (16) and (17), we can eliminate ϕ to obtain the trace and determinant for equilibria,
det ( J ) = ε ( κ ( κ 1 ) / a a ω 2 / κ 2 ) ,
(21)
tr ( J ) = κ ( 1 κ ) / a ε .
(22)

We seek to determine bifurcation curves in ( ω , a )-parameter-space. Our system is planar so that saddle-node and Hopf bifurcations are determined via the conditions det ( J ) = 0 and tr ( J ) = 0 with det ( J ) > 0, respectively. Bogdanov–Takens bifurcations occur at the intersection of Hopf and saddle-node bifurcations. Eliminating ϕ and κ corresponding to equilibria in det ( J ) or tr ( J ) appears not to be feasible since (18) neither can be solved for κ in closed form nor is there a suitable substitution to achieve the elimination. Instead, to determine the desired bifurcation curves, we seek a parameterization of ω and a in terms of one of the equilibrium variables, κ.

1. Saddle-node bifurcation curve
To determine such a parameterization, instead of using the condition det J = 0 for saddle-node bifurcations, we pursue another strategy. Solving the equilibrium condition (18) for
a ( κ ) = ± κ ( κ 1 ) κ 2 ω 2 ,
(23)
and computing
d a d κ = ± ( κ 3 + ω 2 ( 1 2 κ ) ) ( κ 2 ω 2 ) 2 / 3 ,
(24)
we find that d a / d κ = 0 corresponds to a saddle-node condition, resulting in
ω SN ( κ ) = ± κ 3 / 2 2 κ 1 .
(25)

With (25) and (23) [where (25) is inserted for ω], we have obtained the desired parameterization [ ω SN ( κ ) , a SN ( κ )] for the saddle-node bifurcation.

2. Hopf bifurcation curve
Solving tr ( J ) = 0 for a in (22), we find that
a H ( κ ) = κ ( 1 κ ) / ε .
(26)
From (18), we obtain the equilibrium value for
ω H ( κ ) = κ a a 2 κ 2 + 2 κ 1 .
(27)
Hopf curves are, thus, determined via (26) and (27) with κ as parameter, provided that det ( J ) in (21) is positive. Note that all Hopf bifurcations in this system are subcritical so the limit cycles created in the Hopf bifurcations are unstable.
3. Homoclinic bifurcation of (librational) limit cycles

When a librational limit cycle (born in a Hopf bifurcation) collides with a saddle point, it may eventually get destroyed in a homoclinic bifurcation. Homoclinic bifurcations have been determined numerically as follows. The associated homoclinic bifurcation has been determined by inspection of phase portraits while varying parameters, i.e., for a fixed value of a by determining the value of ω at which the limit cycle is destroyed in saddle collision. When limit cycles in this planar system are unstable, Eqs. (15) have been numerically integrated backwards in time while continuing ω. Doing so for several values of a allowed to construct the associated homoclinic bifurcation curve shown in Fig. 3.

4. Bogdanov–Takens points

Curves of saddle-node, Hopf, and homoclinic bifurcations intersect in a codimension 2 Bogdanov–Takens bifurcation point (BT). This point is found by seeking solutions that simultaneously satisfy equilibrium conditions together with det ( J ) = tr ( J ) = 0.

5. Cusp bifurcation points

Another codimension 2 bifurcation point is the cusp bifurcation (CP) where two saddle-node bifurcations meet. Thus, we find that CP points are located at ( ω CP , a CP ) = ( 0 , ± 1 ).

6. Heteroclinic bifurcation and loss of stability of drift cycle

A similar procedure as in Sec. III B 2 c was used to determine bifurcations of the drift cycle, the rotational limit cycle (limit cycle wrapping around the cylinder). This drift cycle is destroyed in a (global) heteroclinic bifurcation (Het, see Fig. 3). However, for a < 0, the drift cycle loses appears to lose its stability (D-loss, see Fig. 3) before a saddle collision may occur; instead, trajectories approach a stable spiral.

3. Stability and bifurcation diagrams

We explain the bifurcation scenarios and associated stability boundaries based on our previous analysis. We first observe that the bifurcation and stability diagram (see Figs. 2 and 3, respectively) exhibit certain symmetries regarding the reflection ω ω. It is easy to check that Eqs. (15), (18), and (19) obey the symmetry (A4), which effectively swaps the labels of the two oscillators and, thus, preserves equilibria, their stabilities, and bifurcations (as seen in Figs. 2 and 3). A further symmetry regarding a a (A5) preserves equilibria, but the stability of the equilibria is not preserved as can be seen by inspecting the Jacobian in (20) [more specifically, saddle-node and cusp bifurcations are preserved, while Hopf bifurcations are not, see Figs. 3, 2(b), and 2(c)]. See also the  Appendix for further details and explanations.

The non-adaptive case ( a = 0) represents the limit of the classical Kuramoto model (see Sec. III A). This case illustrates the perhaps most fundamental bifurcation organizing the structure of the stability diagram on a “global” scale. When intrinsic frequencies are sufficiently similar, ω = ω 1 ω 2 < 1, oscillators lock their frequency, i.e., synchronization is achieved; vice versa, if intrinsic frequencies are too dissimilar, ω > 1, oscillator phases drift apart. Close to ω = 0, oscillators are nearly in-phase (compare with Fig. 2 where a 0); but the phase difference ϕ increases with larger | ω | until the locking breaks apart in a saddle-node bifurcation (black curve) and a drifting state appears. This saddle-node bifurcation separates the stability regions for frequency-locked (L, light gray shaded region) and drifting (D, white region) solutions and can be followed throughout the stability diagram in Fig. 3. Note that this drift solution corresponds to a rotation on the cylinder T × R.

For non-zero adaptivity ( a 0), the structure of the stability diagram in Fig. 3 is organized around two bifurcation points of codimension two, from where additional bifurcation curves emanate: (i) cusp points (CP, black dot) located at ( ω CP , a CP ) = ( 0 , ± 1 ) and (ii) Bogdanov–Takens points (BT, purple dot). Furthermore, several parameter regions of bistability (dark gray shaded regions, yellow regions) appear near the cusp bifurcation points (dark gray shaded regions) and the saddle-node bifurcations leading to drifting solutions (yellow regions).

1. Positive adaptivity (a > 0)

For smaller frequency mismatch (around | ω | < 1) and intermediate adaptivity ( 0 < a < 1), frequency-locked solutions (L, light gray shaded region in Fig. 3) with positive coupling strengths exist, stable nodes [black curve in Fig. 2(b)] closer to in-phase and saddles (black dashed curve) closer to anti-phase configurations. For a 1, two additional saddle-node bifurcation curves [black dots in Fig. 2(a)/black curves in Fig. 3] emerge in a cusp bifurcation (CP, black dot in Fig. 3) located at ( ω CP , a CP ) = ( 0 , 1 ). Traversing these saddle-node bifurcation curves gives birth to an additional stable node and a saddle with negative coupling, see Fig. 2(a). Thus, for a > 1, above the cusp point, between the inner SN bifurcation curves, we have a region of bistability (dark gray shaded region) with two stable locked solutions characterized by near in- and anti-phase configurations, respectively.

For larger frequency mismatch (around | ω | > 1), as we diminish | ω | and traverse the saddle-node bifurcation giving rise to the lock state (L), we first observe a region of bistability (D/L, yellow region in Fig. 3) between the lock state (L) and a drift cycle (D). As the frequency mismatch | ω | is further diminished, the drift cycle, which already was present in the region D for larger | ω | collides with the saddle version of the lock state emanating from the saddle-node bifurcation and is, thus, destroyed in a heteroclinic bifurcation. (labeled Het in Fig. 3; see also Sec. III B 2 c). Note this heteroclinic bifurcation is global and expected to depend on ϵ.

2. Negative adaptivity (a < 0)

The symmetry (A5) swaps near in-phase with anti-phase configurations and changes the stability of equilibria. The overall bifurcation structure is, thus, similar to the one observed for a > 0; however, their structure is more intricate as explained in the following.

First we restrict our attention to smaller frequency mismatch ( | ω | < 1). For 1 < a < 0, the saddle-node bifurcations [black dot in Fig. 2(c)/black curve in Fig. 3] now give birth to a saddle and an unstable node. The unstable node turns into an unstable spiral, which later gains stability in a Hopf bifurcation (HB). The resulting stable spiral [red line in Fig. 2(c)] now, however, corresponds to a near in-phase state with positive coupling (but smaller compared to a > 0). We observe a cusp bifurcation point (CP, black dot in Fig. 3) at ( ω CP , a CP ) = ( 0 , 1 ). Below the cusp point ( a < 1), a Bogdanov–Takens point (BT 2, purple dot in Fig. 3) appears on the associated saddle-node bifurcation [black dot in Fig. 2(d)/black curve in Fig. 3], which gives rise to a subcritical Hopf bifurcation [blue dot in Fig. 2(d)/blue curve in Fig. 3] that stabilizes the second lock state, a stable spiral with small negative coupling and ϕ near ± π / 4 [red line in Fig. 2(d)]. The unstable limit cycles (librations) emanating from the Hopf bifurcation get destroyed in a homoclinic bifurcation (purple curve).

For larger frequency mismatch ( ω > 1), we observe a region of bistability (D/L, yellow region in Fig. 3) between lock (L) and drift (D) states similar to the one observed for positive adaptivity, a > 0. However, there are important differences. First, below the Bogdanov–Takens point BT 1, the frequency-locked solution (L) loses stability in a subcritical Hopf bifurcation (blue curve in Fig. 3). More specifically, as shown in Figs. 2(c) and 2(d), the saddle-node bifurcation (black dot) gives birth to a saddle (black dashed curve) and an unstable node, which first becomes an unstable spiral (red dashed curve) and then becomes a stable spiral (red line) in a Hopf bifurcation (blue dot).

Second, for adaptivity values roughly above the BT 1 point, [the precise location of the origin for this boundary remains unclear (see also detection method in Sec. III B 2)] the drift cycle D is destroyed in a heteroclinic bifurcation; this scenario, however, appears to be different for adaptivity a below BT 1: On the boundary of the D/L (labeled D-loss in Fig. 3; see also Sec. III B 2 f) and L regions (light gray shaded region)—before any heteroclinic bifurcation occurs—the drift (limit) cycle appears to lose its stability. Thus, trajectories emanating from the unstable drift end up spiraling into the stable spiral that arises in the Hopf bifurcation related to BT 1.

To summarize the respective differences regarding stability changes occurring for a < 0 and a > 0, the presence of Bogdanov–Takens bifurcation points BT 1 and BT 2 diminishes the size of the L locking region and the 2L bistability region (dark gray shaded region in Fig. 3). Furthermore, the size of the bistable D/L region with lock and drift is also effectively diminished.

We now allow for arbitrary values of β; in general, this may lead to asymmetric coupling strengths, κ 12 κ 21, corresponding to directed network weights.

1. Stability analysis

We first investigate equilibria of the full three-dimensional system in Eq. (7), which gives rise to the following fixed point conditions:
2 ω = κ 21 sin ( α + ϕ ) κ 12 sin ( α ϕ ) ,
(28a)
κ 12 = 1 + a cos ( β + ϕ ) ,
(28b)
κ 21 = 1 + a cos ( β ϕ ) .
(28c)
Eliminating κ 12 and κ 21, we obtain a condition that only depends on ϕ and guarantees the existence of an equilibrium,
ω = ( cos α + a cos ( α β ) cos ϕ ) sin ϕ .
(29)
Note that Eqs. (28b), (28c), and (29) provide us with a parameterization of an equilibrium curve for ( ω , κ 12 , κ 21 ) in variable ϕ.
To determine the stability of equilibria, we calculate the Jacobian for (7),
J = [ J ϕ ϕ 1 2 sin ( α ϕ ) 1 2 sin ( α + ϕ ) ε a sin ( β + ϕ ) ε 0 ε a sin ( β ϕ ) 0 ε ] ,
where J ϕ ϕ = ( κ 12 cos ( ϕ α ) + κ 21 cos ( ϕ + α ) ) / 2. The Jacobian has the eigenvalues
λ 1 , 2 = μ ± δ ,
(30a)
λ 3 = ε ,
(30b)
where
μ = A 4 ε 2 ,
(31a)
δ = ε 2 ( a B A ) + 1 16 ( A + 2 ε ) 2 ,
(31b)
A := κ 12 cos ( α ϕ ) + κ 21 cos ( α + ϕ ) ,
(31c)
B := cos ( α + β ) cos 2 ϕ cos ( α β ) .
(31d)
Since by definition λ 3 < 0, it suffices to only consider λ 1 , 2.
If an equilibrium point satisfies μ 2 = δ, one eigenvalue becomes zero, and the equilibrium is a saddle-node point. One could use this condition to determine the associated saddle-node curves, but we instead use the following consideration. Since the equilibrium condition (29) only depends on ϕ, the conditions for an equilibrium and for a saddle-node bifurcation are reduced to a problem in a single variable, ϕ. Accordingly, we require in addition to (29) that d ω / d ϕ = 0. Solving the resulting two conditions for ( ω , a ) results in a parameterization of the saddle-node curves in ϕ,
ω SN = cos α sec 2 ϕ sin 3 ϕ ,
(32a)
a SN = cos α cos ϕ sec ( α β ) sec 2 ϕ .
(32b)
Next, we consider equilibria undergoing Hopf bifurcations. If a given equilibrium satisfies μ = 0 , δ < 0, the real part of a complex conjugated pair of eigenvalues vanishes and the equilibrium is a Hopf point. We find a parameterization for the Hopf curves as follows. First, we require that the equilibrium condition in (29) is satisfied. Second, after substituting (28b) and (28c) into (31c), we require that μ = 0. Solving the two resulting conditions for ( ω , a ), we obtain
ω H = ( ε cos α cos β cos ϕ + ( cos α + ε cos ϕ ) sin α sin β ) sin ϕ sin α sin β sin 2 ϕ cos α cos β cos 2 ϕ ,
(33a)
a H = ε + cos α cos ϕ sin α sin β sin 2 ϕ cos α cos β cos 2 ϕ .
(33b)
Now ( ω H , a H ) with parameter ϕ delineates a Hopf curve provided that δ < 0 in (31b).

Finally, if μ = δ = 0, two eigenvalues are zero and the equilibrium is a Bogdanov–Takens point. Thus, simultaneously solving μ = δ = 0 together with the fixed point conditions (28), given a value of ε we are able to obtain ( ϕ , κ 12 , κ 21 , a , ω ) at a Bogdanov–Takens point.

All other bifurcations were detected using numerical continuation (MatCont46), except for homoclinic and heteroclinic bifurcation boundaries, and loss of stability of the drift cycle (see also Sec. III B 2 c).

2. Dynamics for β=π/2

Substituting β = π / 2 into Eqs. (9), we have
d Δ d t = ε ( 2 a sin ϕ + Δ ) ,
(34a)
d Σ d t = ε ( 2 Σ ) ,
(34b)
so that Σ 2 as t and the system is attracted to a two-dimensional subspace on which the dynamics are given by
d ϕ d t = ω + sin α 2 Δ cos ϕ cos α sin ϕ ,
(35a)
d Δ d t = ε ( 2 a sin ϕ + Δ ) .
(35b)
Note that in the asymptotic time limit t the coupling is not anti-symmetric, i.e., we have κ 12 = κ 21; this would be the case if instead we had chosen a 0 = 0 in our model Eq. (4). Instead, we have κ 12 + κ 21 2.
1. Phase-lag α = 0 and α = π
We first discuss the simplest case with α = 0. The asymptotic dynamics in
d ϕ d t = ω sin ϕ ,
(36a)
d Δ d t = ε ( 2 a sin ϕ + Δ ) .
(36b)
While the variable ϕ drives Δ, its dynamics is independent of Δ. Thus, we effectively observe the one-dimensional dynamics known for the Kuramoto model (11) discussed in Sec. III A.
For ω < 1, we observe locked states with ω = sin ϕ. The resulting equilibria are
ϕ = { π arcsin ω , arcsin ω ,
(37a)
κ 12 = 1 a ω ,
(37b)
κ 21 = 1 + a ω .
(37c)
where the latter near in-phase state is stable. The equilibrium condition, ω = sin ϕ, informs us that a saddle-node bifurcation occurs for ω = ± 1, regardless of the value of a. The associated stability boundaries are, therefore, the straight black curves shown in Fig. 4(a). Thus, contrary to other parameter choices considered in this study, the dynamics and bifurcations do not increase in complexity as a is varied: we are dealing with a special limiting case. The equilibria are the same for α = π, however, with reversed stability.
2. Arbitrary phase-lag (α ≠ 0, π)

It is easily checked that (34a) has the parameter symmetry ( a , α ) ( a , α ). Since we were allowed to limit α ( π / 2 , π / 2 ] in the original governing Eq. (7), we may restrict α further to either [ π / 2 , 0 ] or [ 0 , π / 2 ] while observing this symmetry.

To keep the analysis manageable, we consider only α = π / 10. The resulting stability diagram is shown in Fig. 4(b). Even though this is just a small deviation, the bifurcation landscape differs drastically from the case where α = 0. As was the case for α = β = 0 (Fig. 3), the stability diagram is symmetric regarding ω ω [preserving equilibria and their stabilities and bifurcations, see (A4)]. Moreover, equilibria and the SN and cusp bifurcations are symmetric about a a [preserving equilibria and the SN and cusp bifurcations, see (A5)]. See the  Appendix for further details and explanations.

3. Dynamics for β=π/4

Next, we consider the case where β = π / 4. The stability diagram for α = 0 (not shown) is qualitatively identical to the one obtained for α = β = 0 (Fig. 3). Therefore, we instead consider small deviations from α = 0, i.e., α = ± π / 10.

1. Phase-lag α = π/10

First, we consider the stability diagram for α = π / 10, see Fig. 5. There are three distinct Bogdanov–Takens points [BT 1 and BT 2 in the main plot of Fig. 5 and (two) BT 3 in the inset], which organize the bifurcation structure. The Hopf curves (blue curves) emanating from all three Bogdanov–Takens points (BT 1, BT 2, BT 3) are subcritical. These Hopf curves are adjacent to SN curves (black curves) giving birth to saddles and unstable nodes. The unstable nodes turn into unstable spirals, and, when undergoing the Hopf bifurcations, into stable spirals. The homoclinic curves (purple curves) adjacent to the SN curves destroy the unstable (libration) limit cycles created in the subcritical Hopf bifurcations.

2. Phase-lag α = −π/10

Second, we consider the stability diagram for α = π / 10 in Fig. 6. The resulting dynamics become far more involved when compared to the previously considered cases. The Hopf curve (blue curve) that emanates from the BT 1 point [purple dot in Fig. 6(a)] contains a Generalized Hopf (GH) point located at ( ω GH, a GH) [blue dot in Fig. 6(a)], which separates the Hopf curve (blue curve) into a supercritical segment (below GH) and a subcritical segment (above GH).

1. Regime of subcritical Hopf bifurcations (a > aGH)

The subcritical Hopf curve [blue curve in Fig. 6(a)] produces unstable (libration) limit cycles, which may undergo a saddle-node-of limit cycles (SNLC) bifurcation (black dashed curve) before they die in a homoclinic bifurcation (purple curve). The homoclinic bifurcation curve, which originates in the point BT 1 (purple dot), was found via numerical continuation using MatCont.46 Below the GH point (blue dot), there are two branches of SNLC bifurcations which meet in a cusp of SNLCs (CPC, green dot). Thus, between the GH and the CPC points, limit cycles produced in the Hopf bifurcation (blue curve) undergo two SNLC bifurcations before dying in the homoclinic bifurcation (purple curve). One of the SNLC branches (black dashed curve) meets the homoclinic bifurcation curve (purple curve) in a point which we call SLH (red dot).

2. Regime of supercritical Hopf bifurcations (a < aGH)

The supercritical Hopf bifurcations [blue curve in Fig. 6(a)] below the GH (blue dot) produce stable (libration) limit cycles (LB), which can only exist on the same side of the supercritical Hopf curve as the unstable spirals, i.e., in the blue region between the HB curve (blue) and the HC curve (purple) that destroys the limit cycles. However, the stable limit cycle can be destroyed/destabilized before reaching the HC bifurcation (purple curve) in an intriguing bifurcation scenario. We let a = 1.9425 and follow the green line in the diagram [Fig. 6(a)] while varying ω and show example trajectories for chosen values of ω in Fig. 7.

FIG. 4.

Stability diagrams for β = π / 2 where ϵ = 0.2. ( ω , a ) parameter regions are shaded according to which stable states they contain Drift (D) and/or locked states (L).

FIG. 4.

Stability diagrams for β = π / 2 where ϵ = 0.2. ( ω , a ) parameter regions are shaded according to which stable states they contain Drift (D) and/or locked states (L).

Close modal
FIG. 5.

Stability diagram for α = π / 10 , β = π / 4 with ϵ = 0.2. ( ω , a ) parameter regions are shaded according to which stable states they contain Drift (D) and/or locked states (L).

FIG. 5.

Stability diagram for α = π / 10 , β = π / 4 with ϵ = 0.2. ( ω , a ) parameter regions are shaded according to which stable states they contain Drift (D) and/or locked states (L).

Close modal
FIG. 6.

Dynamics for α = π / 10 , β = π / 4 and ϵ = 0.2. (a) Stability diagram. ( ω , a ) parameter regions are shaded according to which stable states they contain Drift (D, pure rotation around cylinder) and/or locked states (L), librations (LB, periodic or chaotic), mixed (libration/rotation) oscillation (MO, periodic or chaotic). (b) Bifurcation diagram for α = π / 10 , β = π / 4 , a = 1.9425 [value of a is marked as green line in (a)]. The diagram shows stable nodes (black curves), saddles, unstable nodes (both black dashed curves), stable spirals (red lines), unstable spirals (red dashed curves), stable oscillations (local extrema marked in blue), and unstable limit cycles (local extrema marked as blue dashed curves). (c) Magnification of (b) [see green rectangle in (b)]. The oscillations (LB/MO) alternate between periodic and chaotic as ω increases, and for ω > 0.782 99 the oscillations are exclusively pure drift states (D). (d): magnification of (c) [see green rectangle in (c)]. The first PD bifurcation (left green dashed line) produces an unstable limit cycle (blue dashed curve), which is destroyed in a homoclinic bifurcation (HC), and a stable, period-2 limit cycle (LB) which undergoes a second PD bifurcation (right green dashed line) and then a cascade of PD bifurcations (not marked, for readability), leading to stable chaos. An MO state co-existing with the limit cycle is marked in magenta for distinguishability. At ω 0.6955, the stable chaotic attractor turns into a stable periodic oscillation. The maximum Lyapunov exponent λ max of the stable oscillatory states is shown beneath. For further explanations, see the text. For phase portraits of the oscillations, see Fig. 7.

FIG. 6.

Dynamics for α = π / 10 , β = π / 4 and ϵ = 0.2. (a) Stability diagram. ( ω , a ) parameter regions are shaded according to which stable states they contain Drift (D, pure rotation around cylinder) and/or locked states (L), librations (LB, periodic or chaotic), mixed (libration/rotation) oscillation (MO, periodic or chaotic). (b) Bifurcation diagram for α = π / 10 , β = π / 4 , a = 1.9425 [value of a is marked as green line in (a)]. The diagram shows stable nodes (black curves), saddles, unstable nodes (both black dashed curves), stable spirals (red lines), unstable spirals (red dashed curves), stable oscillations (local extrema marked in blue), and unstable limit cycles (local extrema marked as blue dashed curves). (c) Magnification of (b) [see green rectangle in (b)]. The oscillations (LB/MO) alternate between periodic and chaotic as ω increases, and for ω > 0.782 99 the oscillations are exclusively pure drift states (D). (d): magnification of (c) [see green rectangle in (c)]. The first PD bifurcation (left green dashed line) produces an unstable limit cycle (blue dashed curve), which is destroyed in a homoclinic bifurcation (HC), and a stable, period-2 limit cycle (LB) which undergoes a second PD bifurcation (right green dashed line) and then a cascade of PD bifurcations (not marked, for readability), leading to stable chaos. An MO state co-existing with the limit cycle is marked in magenta for distinguishability. At ω 0.6955, the stable chaotic attractor turns into a stable periodic oscillation. The maximum Lyapunov exponent λ max of the stable oscillatory states is shown beneath. For further explanations, see the text. For phase portraits of the oscillations, see Fig. 7.

Close modal
FIG. 7.

Regime of supercritical Hopf bifurcations ( a < a GH): phase portraits of stable oscillations for α = π / 10 , β = π / 4 , a = 1.9425, ϵ = 0.2 and varying ω [shown as green lines in Figs. 6(b)6(d)]. (a) and (b) Period-1 (a) or period-2 (b) libration cycle (LB, blue) co-exists with mixed rotational/librational drift cycle (MO, magenta). (c) Chaotic libration (LB, blue). (d) Chaotic mixed oscillation (MO), consisting of rotations and librations. (e) Periodic MO with 1 rotation and 4 librations (blue). (f) Periodic MO with 2 rotations and 1 libration. (g) Periodic MO with 3 rotations and 1 libration. (h) Pure rotation drift cycle (D, blue). See the text for further details.

FIG. 7.

Regime of supercritical Hopf bifurcations ( a < a GH): phase portraits of stable oscillations for α = π / 10 , β = π / 4 , a = 1.9425, ϵ = 0.2 and varying ω [shown as green lines in Figs. 6(b)6(d)]. (a) and (b) Period-1 (a) or period-2 (b) libration cycle (LB, blue) co-exists with mixed rotational/librational drift cycle (MO, magenta). (c) Chaotic libration (LB, blue). (d) Chaotic mixed oscillation (MO), consisting of rotations and librations. (e) Periodic MO with 1 rotation and 4 librations (blue). (f) Periodic MO with 2 rotations and 1 libration. (g) Periodic MO with 3 rotations and 1 libration. (h) Pure rotation drift cycle (D, blue). See the text for further details.

Close modal

After the supercritical Hopf bifurcation [right blue dot in Fig. 6(b)], i.e., for increasing ω, the stable limit cycle LB [blue in Figs. 6(b), 6(d), and 7(a)] undergoes a Period-Doubling (PD) bifurcation [left green dashed line in Fig. 6(d)]. The resulting stable limit cycle [blue in Figs. 6(d) and 7(b)] has period 2. The original, period-1 cycle has become unstable [blue dashed curve in Fig. 6(d)] in the PD [left green dashed line in Fig. 6(d)]; the associated branch is ultimately destroyed in the homoclinic bifurcation [purple dot in Figs. 6(b) and 6(c)]. Remarkably, note that the libration cycle LB co-exists with a (rotational) drift limit cycle in the region of period-1 and period-2 cycles, see magenta curves in Fig. 6(d) and magenta trajectories in Figs. 7(a) and 7(b). Unlike the previously described D state, these drift cycles consist of both rotations and librations, which is why refer to them as mixed oscillations (MOs). While the precise nature of the bifurcation mechanism giving rise to its birth on the left for smaller ω remains open, inspection of trajectories for larger ω on the right strongly suggests that this mixed oscillation is destroyed in a collision with the unstable period-1 version of the libration cycle [blue dashed curve in Fig. 6(d)].

As we increase ω further, the stable period-2 cycle undergoes a period-doubling cascade (PDs are not explicitly marked to improve readability) resulting in a chaotic attractor with ϕ bounded, i.e., it can (still) be seen as a libration LB [blue in Figs. 6(c), 6(d), and 7(c)]. To characterize this chaotic motion, we numerically estimated the maximum Lyapunov exponent λ max, and we found that λ max > 0 (close to zero) throughout the (non-)chaotic ranges of ω [Figs. 6(c) and 6(d)]. (Fluctuations in λ max while varying ω are due to the fact that phase points of the original and the perturbed trajectory can lie anywhere on the chaotic attractor by the simulation time λ max is computed.)

This chaotic attractor can be seen as a chaotic libration, since ϕ does not rotate around the cylinder T × R [Fig. 7(c)]. Remarkably, as ω is further increased, the chaotic attractor at some point also includes oscillations characterized as rotations, where ϕ revolves around the cylinder [Fig. 7(d)] (in a fashion that is reminiscent of phase slips). More precisely, the trajectory librates a number of times, until it escapes and becomes a rotation with drift character, while, however, remaining chaotic in nature. Presumably, this behavior is akin to a mixed mode oscillation seen characteristic of slow-fast systems, which is a reminiscence from low ε (recall that we use a fairly large ε = 0.2 here). Eventually, as ω is increased even further, the cycle ceases to display chaoticity and becomes a periodic oscillation including both librations (of a certain period) and (drift) rotations [Fig. 7(e)], which we therefore refer to as mixed oscillations (MOs; not to be confused with mixed-mode oscillations47). Ultimately, as we further increase ω, the librations disappear and the cycle is a purely rotational drift cycle [Fig. 7(h)]. This is remarkable, since—as already noted further above—a stable mixed oscillation cycle was already present for smaller ω [magenta curves in Figs. 6(c), 6(d), 7(a), and 7(b)], which was destroyed in a bifurcation scenario different from the scenario seen here. Thus, these (rotational/librational) drift cycles just described are indeed distinct, as they seem to emerge from the libration cycles originally born in the supercritical Hopf bifurcation. To summarize, stable oscillations encountered are all either non-drifting libration cycles (LB, periodic or chaotic), purely rotational drift states (D, periodic), or mixed rotational/librational drift oscillations (MO, periodic or chaotic).

Finally, note that states containing both librations and rotations exist, for a = 1.9425, in the range ω [ 0.694 56 , 0.782 99 ) [Figs. 6(c) and 7(d)7(g)], where subintervals hosting periodic cycles alternate with subintervals hosting (aperiodic) chaotic motion. This was found via quasi-continuation of the stable oscillatory states in the interval ω [ 0.693 , 1.25 ]. Interestingly, the respective number of rotations and librations per period of the oscillation is different for each periodic subinterval [see Figs. 7(e)7(g)]. Finally, based on quasi-continuation of the drift states in the direction of ω 0 for several values of a, we find that stable purely rotational drift states (D) exist only in the white regions and yellow regions in Fig. 6(a). While the blue region may contain all aforementioned types of LB and MO cycles, none of these types emerges everywhere inside of the blue region; for readability, we do not distinguish existence regions for each type individually.

In Sec. III, we have studied in detail the behavior of the system (5) for N = 2 oscillators. The question abounds whether some of the behavior observed for only N = 2 oscillators carries over to a larger version of the system with N = 50 oscillators. In order to address this question, we restrict our attention to the simplest case with parameters α = β = 0 , ε = 0.2, (see Figs. 2 and 3 for the results obtained for N = 2). For two oscillators, we could tune the frequency mismatch, ω; for the many oscillator systems, we instead let the intrinsic frequencies ω l follow a normal distribution with zero mean and standard deviation σ. Furthermore, instead of drawing these frequencies randomly, we constructed them via an equidistant set of points mapped to the interval [ 3 σ , 3 σ ] via the inverse error function. This symmetrization eliminates potential random effects that due to finite size effects might otherwise obscure the dynamic effects we are interested in, such as to render comparison to the case of N = 2 more immediate.

We study the dynamic asymptotic behavior of this system on a grid in ( σ , a ) parameter space, the results of which are summarized in Fig. 8(a). We may distinguish four states: the (frequency-)locked state (L) and the antipodal state (AP), two frequency-locked states where all d ϕ l / d t , d κ l m / d t tend to zero after some transient; a drift (D) state with all oscillator phases drifting a part; and the partial synchrony state (PS) which can be seen as a mixture of both L and D states, in the sense that the N oscillators are split into a frequency-locked and a drifting subset. We describe these states in the order of appearance as we increase σ for a = 1.7 as shown in Fig. 8.

FIG. 8.

Summary of simulation results for N = 50 oscillators with α = β = 0 , ε = 0.2. (a) Stability diagram for varying a and σ summarizes the four post-transient final states (as observation for one realization only): antipodal (AP, triangle), frequency-locked (L, black circles), partially synchronous (PS, gray circles), drift (D, empty circles). Panels (b)–(d) display the behavior observed for three distinct parameter values [marked by vertical arrows in (a)], i.e., AP for a = 1.7 , σ = 0.2 (b), L for a = 1.7 , σ = 0.4 (c), and PS a = 1.7 , σ = 0.6 (d), respectively. The top row in the panels displays phases (blue/red dots for locked/drifting oscillators), ϕ l, as well as dynamic frequencies (black dots), d ϕ l d t in Eq. (5), at the final simulation time T. The bottom row in panels (b) and (c) features the κ matrices at final time T = 1500, where white and black pixels correspond to 1 a = 0.7 and 1 + a = 2.7, respectively (both matrices have zero diagonal shown in gray). Note that the black/white regions are not uniformly valued, but heterogeneous values (see the text). The bottom row in panel (d) presents the time evolution of κ l m ( t ) during 0 t T = 200 (transients have subsided for t 50).

FIG. 8.

Summary of simulation results for N = 50 oscillators with α = β = 0 , ε = 0.2. (a) Stability diagram for varying a and σ summarizes the four post-transient final states (as observation for one realization only): antipodal (AP, triangle), frequency-locked (L, black circles), partially synchronous (PS, gray circles), drift (D, empty circles). Panels (b)–(d) display the behavior observed for three distinct parameter values [marked by vertical arrows in (a)], i.e., AP for a = 1.7 , σ = 0.2 (b), L for a = 1.7 , σ = 0.4 (c), and PS a = 1.7 , σ = 0.6 (d), respectively. The top row in the panels displays phases (blue/red dots for locked/drifting oscillators), ϕ l, as well as dynamic frequencies (black dots), d ϕ l d t in Eq. (5), at the final simulation time T. The bottom row in panels (b) and (c) features the κ matrices at final time T = 1500, where white and black pixels correspond to 1 a = 0.7 and 1 + a = 2.7, respectively (both matrices have zero diagonal shown in gray). Note that the black/white regions are not uniformly valued, but heterogeneous values (see the text). The bottom row in panel (d) presents the time evolution of κ l m ( t ) during 0 t T = 200 (transients have subsided for t 50).

Close modal

For the antipodal/anti-phase state (AP) [Fig. 8(b)], the phases of oscillators split into two groups inside which phases differ only very little from one another, but in between the groups phases differ by about π, i.e., the two oscillator groups are in an anti-phase (antipodal) configuration with respect to one another on the phase circle; d ϕ l / d t , d κ l m / d t tend to zero asymptotically in time. This special configuration of the phases also impacts the coupling on the network: κ l m that represent couplings between oscillator pairs residing in the same population (distinct populations) assume values close to the maximum possible value 1 + a (the minimum possible value 1 a), see black (white) κ l m values in Fig. 8(b). Note that these values for κ l m are heterogeneously distributed near these values as they accommodate for the non-identical phase configuration (the figure compresses the variation due to the chosen gray scale). Note that the antipodal state occurs when the distribution of ω l is sufficiently narrow and the adaptivity a > 0 is sufficiently strong. This state may be interpreted as a corollary to the 2L co-existence observed for N = 2 , α = β = 0 , a 0 , | ω | 1 [see Fig. 2(a)], where the two oscillators either occupy a phase configuration with 0 or π difference; here, on the other hand, oscillators form groups of varying size that mutually adhere to either of the two configurations. Note that a similar correspondence applies to the coupling strengths both at the upper and the lower end of the applicable κ range.

For the frequency-locked state L [Fig. 8(c)], all phases are close to one another. This configuration implies that the κ l m are close to 1 + a (and bounded by that value from above). Note that individual coupling strengths κ l m for l and m are not identical but heterogeneously distributed: the coupling strengths have adapted to accommodate for the non-identical phase configuration. In Figs. 8(b) and 8(c), the gray scale for the values of κ l m ranges from 1 a to 1 + a so that the distributed values cannot easily be distinguished. Topologically, this essentially corresponds to an all-to-all network. This state is, thus, interpreted as the corollary of the frequency-locked state L for N = 2, see Fig. 2(b) for small | ω |.

For the partially synchronous state (PS) [in Fig. 8(d)], oscillators split into two groups where one is locked (blue dots for phases), while the other contains drifters (red dots for phases). Correspondingly, coupling strengths κ l m where l and/or m can be associated with drifting oscillators and are, thus, oscillatory [red curves in Fig. 8(d)]. By contrast, the locked oscillators [blue dots in Fig. 8(d)] display vanishingly small dynamic frequencies; thus, the κ l m for which neither l nor m are drifting oscillate close to 1 + a and with very small amplitude [blue time traces in Fig. 8(d)]. We note that the correspondence to states observed for N = 2 is less obvious; however, an extension of the D/L bistable region for N = 2 oscillators to larger systems could give space to accommodate for configurations where a group of oscillators populate a lock state while the remaining oscillators adhere to drifting oscillations.

For the drifting state (D), oscillator phases are incoherent and drift apart from each other in unlocked motion while the coupling strengths κ l m appear to attain a (quasi-)stationary distribution centered around 1 (not shown). We interpret this state as the corollary to the D state for N = 2.

Just as for the N = 2 state, the locked state L occurs for both positive and negative a. For N = 50 and a < 0 with σ small, we observe locked states (corresponding to stable foci for N = 2) whose κ l m are mostly close to the minimal value 1 a, in analogy to the case of N = 2 oscillators [see Fig. 2(c)].

For the N = 2 system, we found that (co)existence of one (L) or two (2L) stable locked state(s) is guaranteed for | a | sufficiently large. Interestingly, Fig. 8(a) suggests that there is a stable locked state for a > 0 sufficiently large; yet, for a < 0 it appears as if the region of stable locked states is diminished in favor of partially synchronous, and, even drift states (incoherence). Similarly, we observed for N = 2 and a < 0 that the locking region is diminished when compared to a > 0, which is due to a subcritical Hopf bifurcation that destabilizing the frequency-locked branch (albeit not as much as strongly as for N = 50)—we hypothesize that this effect is more pronounced for larger oscillator numbers (and ensuing bifurcations). It would be interesting to investigate this aspect in future research.

The model in Eq. (5) is a Kuramoto–Sakaguchi model with phase-lag α and coupling strengths that adapt according to a learning rule with adaptivity strength a and adaptation shift β. The parameter a allows to tune away from the limit of the non-adaptive classical Kuramoto model with stationary coupling ( a = 0) and systematically investigate the effect of adaptivity on the collective dynamics.

Our bifurcation analysis concerns the case of N = 2 oscillators, for which we—chiefly—may distinguish three paradigms of adaptivity: (i) the non-adaptive limit with stationary coupling ( a = 0); (ii) symmetric adaptation ( β = 0 or β = π); and (iii) asymmetric adaptation (arbitrary β). The non-adaptive limit ( a = 0) trivially reduces to the classical Kuramoto model, where, for frequency locking to occur, the frequency mismatch must be smaller than the coupling strength between oscillators.

Considering the case of symmetric coupling ( β = 0), we first observe that deviating from the non-adaptive limit ( a = 0) with non-zero adaptivity ( a 0) leads to an overall larger locking region (L), i.e., a larger frequency mismatch is required to break locking. It is interesting to note that in the context of forced oscillations, the mode-locking region corresponds to an Arnol’d tongue; indeed, the strength of adaptivity could be related to a forcing strength acting on one of the oscillators, at least in the limit of slow adaptation ( ε 0). Furthermore, smaller frequency mismatch | ω | and sufficiently large adaptivity a allow for bistable regions 2L where two frequency locked modes exist, i.e., anti-phase (antipodal) configurations co-exist in addition to in-phase configurations—indeed, considering stationary coupling in, e.g., Eq. (15) effectively results in a bi-harmonic coupling function which may lead to such bi-stability.48 For larger | ω | and a 0, a further bistability region D/L appears where locked states co-exist with drift cycles (rotations). While this general picture prevails for both positive and negative adaptivity, the situation is slightly more complicated in the region with a < 0, where additional bifurcations diminish the width of the various (bi-)stable regions (L, 2L  and D/L). Drift cycles correspond to rotational limit cycles around the cylindrical phase space, whereas Hopf bifurcations give rise to librational limit cycles; however, in the case of symmetric coupling ( β = 0 , π), these cycles remain always unstable.

We found that similar bifurcation scenarios occur for networks with asymmetric coupling for β 0 and a > 0; but, even more intriguing dynamics are possible for a < 0 and β = π / 4 , α = π / 10. In short, a generalized Hopf bifurcation (GH) may be present, around which a complex structure of bifurcations is organized [including sub-/supercritical Hopf, saddle-node of limit cycles (SNLC), cusp of cycles (CPC), a homoclinic and SLH bifurcations, similar to a scenario seen for a system of Theta neurons, see Ref. 49]; in particular, for a < 0 the GH may give rise to a supercritical Hopf bifurcation leading to stable librational limit cycles. These libration cycles may undergo a further period-doubling cascade to chaos. Moreover, remarkably, period-1 and period-2 libration cycles may co-exist with a rotational drift cycle characterized by a non-trivial winding number; this rotational cycle appears to be destroyed in a collision with the unstable period-1 cycle generated by the first period-doubling bifurcation, before the stable libration cycle becomes chaotic. However, the period-doubling cascade of these librations features a surprise. As ω increases, after the librations have become chaotic, parameter regions exist where the librations are altered into a mixed oscillation, which combines features of both libration and rotation. The winding number of the libration gradually changes, and mixed oscillations alternate between chaotic and periodic, until for large ω one finds a regular period-1 rotational drift cycle.

One may expect that the analysis of the N = 2 oscillator system might be able to capture certain aspects of the dynamic behavior seen in larger systems. Our simulations for adaptive networks with symmetric coupling ( β = 0) with N = 50 revealed that such a correspondence exists, albeit with some limitations. First, considering a > 0, we find a relatively good correspondence between antipodal/anti-phase (AP), locked (L), and drifting/incoherent states (D) seen for N = 50 with the bistable locking region (2L), the locking region (L), and drifting region (D) seen for N = 2; the partially synchronous (PS) state for N = 50 may be related to the D/L region for N = 2. Second, increasing strength of adaptivity | a | leads to a widening of the locking regions for N = 2; this effect appears to be present for N = 50 only when a > 0, but not for a < 0.

We point out similarities between the system studied here and systems subject to previous research. These systems are contained as special cases in our model (5). The relationships are more obvious when considering the unscaled version of our adaptive model (4): Letting a 0 = 0 with α = β = 0, we recover the model Eqs. (3) and (4) studied by Seliger et al., Ref. 8. Similarly, if we set a 0 = 0 with ω 1 = = ω N = 0 , a 1 = 1, we recover the model by Berner et al.10 Our adaptation rule generalizes the model in Ref. 8 by inclusion of an adaptation shift β 0 (see Ref. 10); but in contrast to Ref. 10, we introduced a nontrivial adaptation offset a 0 and strength of adaptation a 1. This choice allowed us to systematically deviate from the classical Kuramoto–Sakaguchi model with stationary coupling, and thereby, to address the question of how varying levels of adaptivity impacts the dynamics of the network. Furthermore, unlike Ref. 10, we allow for non-identical frequencies ( ω 0 and σ 0, respectively) which breaks symmetries in the system—thus, this enables us to study the width of the locking region. We note that the choice a 0 = 1 does not allow us to retrieve the model in Ref. 8 since a 0 0 constitutes a singular limit.

Nevertheless, some of the basic dynamic behavior, such as the presence of frequency locking or antipodal (antiphasic) cluster states seen in Ref. 10, naturally carry over to our system as long as intrinsic frequencies are kept relatively close to each other. Vice versa, breaking the symmetry via disordered frequencies ( ω 0) opens the door for bifurcations resulting in novel states, such as small amplitude librations and the associated period-doubling cascade reported here. A recent study50 carries out an analysis of our system with a 0 = 0 and a 0 = 1 in the slowly adapting limit ( ε 0); however, this analysis separates the dynamics of the coupling strengths into a planar problem, and thus, the resulting dynamics cannot display the transition to chaos or the mixed oscillations containing both librations and rotations that we found here. Finally, we note that Kasatkin and Nekorkin51 studied the same oscillator system with N = 2 as in this study, except that frequencies are identical ( ω = 0); in this case, these authors also observed chaotic dynamics.

Our simulations for N = 50 revealed that oscillators may split into two groups, which experience distinct strong or weak coupling, depending on the specific phase configuration that the oscillators assume. This phenomenon is particularly prominent for the states we called AP and PS. We note that several studies addressed the effect of non-uniform, but time-independent coupling on oscillator networks. In order to understand the occurrence of chimera states, past research performed bifurcation analysis for N = 4 oscillators52 and for larger N, or in the continuum limit ( N ).53–57 However, drawing analogies to our model is not straightforward. First, the configurations seen for the AP or PS state in Fig. 8 display groups with different oscillator numbers, but aforementioned studies consider oscillator groups of equal size. Second, the adaptive capability of the systems renders the dynamics intrinsically more complicated. For instance, the adaptive rule in Eq. (5b) may give rise to bi-harmonic coupling function (depending on the parameter choice) and, thus, to bi-stability of configurations, see Fig. 2, a feature that is due to higher harmonic interaction.48 [This becomes evident, e.g., by assuming a stationary solution for (15).] Moreover, the adaptivity may stabilize certain non-uniform coupling configuration, while models without adaption simply enforce the non-uniform coupling. This is a fundamental difference—as we have seen, adaptive dynamics generates nontrivial dynamic behavior via bifurcations that are absent in the non-adaptive limit (Sec. III A). There exist other models with bimodal frequency distributions where the oscillator population may effectively split into two groups adhering to distinct dynamic behavior,45,58–61 and even models that combine bimodally distributed frequencies with non-uniform coupling strength;62 however, how to correlate the observed states with the behavior observed for these systems seems less obvious.

We plan to further study adaptive Kuramoto–Sakaguchi oscillator networks with large oscillator numbers. To this end, a variety of questions are interesting to pursue. For instance, it will be interesting to see if similar complex dynamics, such as the period-doubling cascade seen for N = 2 with asymmetric coupling ( β 0 , π; see Sec. III C), persists for larger systems. Moreover, developing a (finite) mean-field theory as well as a continuum limit description for Eq. (5), e.g., using graphons/graphops and other methods,33,63–67 would be very useful to study the dynamics of Eq. (5), but remains a formidable challenge. Finally, it would be interesting to study and relate adaptive oscillator networks to more realistic settings such as those that occur in a biological, neuro-scientific, or experimental context.21,25,43 These represent formidable open problems for future research.

We thank M. P. Sørensen, C. Bick, S. Yanchuk, and R. Cestnik for useful discussions and comments throughout this study. We acknowledge the DTU International Graduate School for support via the EU-COFUND project “Synchronization in Co-Evolutionary Network Dynamics (SEND).”

The authors have no conflicts to disclose.

Benjamin Jüttner: Conceptualization (equal); Formal analysis (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal). Erik A. Martens: Conceptualization (equal); Formal analysis (equal); Methodology (equal); Supervision (lead); Writing – original draft (lead); Writing – review & editing (equal).

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

Equation (7) exhibits various symmetries that we survey here. The following symmetries allow us to reduce the α , β parameter range:
( a , β ) ( a , β + π ) ,
(A1)
( a , α , ϕ ) ( a , α + π , ϕ + π ) ,
(A2)
( α , β , κ 12 , κ 21 ) ( α , β , κ 21 , κ 12 ) .
(A3)
Moreover, the following symmetries about the ω , a axes are visible in all stability diagrams:
( ϕ , ω , κ 12 , κ 21 ) ( ϕ , ω , κ 21 , κ 12 ) .
(A4)
As mentioned in the text, this symmetry defined leaves equilibria of Eq. (7), their types, and their stabilities unchanged. Similarly, the symmetry,
( a , ϕ , κ 12 , κ 21 ) ( a , ϕ + π , κ 21 , κ 12 ) ,
(A5)
preserves equilibria of Eq. (7), but the stabilities and bifurcations of these equilibria require further discussion. For this, we look at the eigenvalues of the Jacobian specified in (30). Clearly, the SN condition 0 μ 2 = δ is equivalent to
0 = a B A .
(A6)
If condition (A6) is fulfilled, then it is still fulfilled after applying the symmetry, since the symmetry leaves B unchanged and only inverts the sign of a and A. Looking at Hopf bifurcations, μ = 0 in (31a) is a necessary condition for a Hopf point. If μ = 0 is fulfilled, it is clearly not fulfilled after the symmetry is applied, since the symmetry inverts the sign of A.

In summary, symmetry (A4) preserves equilibria, their stabilities, as well as all bifurcations since we are merely renaming the oscillators; the symmetry (A5) preserves SN and cusp bifurcations, but not Hopf bifurcations. Note that these two symmetries are evident in all stability diagrams shown in this article, i.e., Figs. 3–6(a).

Finally, the following symmetries are not exploited in the text but mentioned here for the sake of completeness:
( α , β , ϕ ) ( α + π , β + π , ϕ + π ) ,
(A7)
( α , ω , ε , t ) ( α + π , ω , ε , t ) ,
(A8)
( a , ω , ε , ϕ , t ) ( a , ω , ε , ϕ + π , t ) ,
(A9)
( α , ω , ε , t ) ( α + π , ω , ε , t ) .
(A10)
1.
S. H.
Strogatz
,
Sync: The Emerging Science of Spontaneous Order
(
Hyperion
,
New York
,
2003
).
2.
A.
Pikovsky
,
M.
Rosenblum
, and
J.
Kurths
,
Synchronization: A Universal Concept in Nonlinear Sciences
(
Cambridge University Press
,
New York
,
2001
).
3.
W.
Singer
and
C. M.
Gray
, “
Visual feature integration and the temporal correlation hypothesis
,”
Ann. Rev. Neurosci.
18
,
555
586
(
1995
).
4.
P. J.
Uhlhaas
and
W.
Singer
, “
Neural synchrony in brain disorders: Relevance for cognitive dysfunctions and pathophysiology
,”
Neuron
52
(
1
),
155
68
(
2006
).
5.
M.
Rohden
,
A.
Sorge
,
M.
Timme
, and
D.
Witthaut
, “
Self-organized synchronization in decentralized power grids
,”
Phys. Rev. Lett.
109
(
6
),
064101
(
2012
).
6.
J.
Acebrón
,
L.
Bonilla
,
C.
Pérez Vicente
,
F.
Ritort
, and
R.
Spigler
, “
The Kuramoto model: A simple paradigm for synchronization phenomena
,”
Rev. Mod. Phys.
77
(
1
),
137
185
(
2005
).
7.
C.
Bick
,
C.
Laing
,
M.
Goodfellow
, and
E. A.
Martens
, “
Understanding the dynamics of biological and neural oscillator networks through exact mean-field reductions: A review
,”
J. Math. Neurosci.
10
),
9
(
2020
).
8.
P.
Seliger
,
S. C.
Young
, and
L. S.
Tsimring
, “
Plasticity and learning in a network of coupled phase oscillators
,”
Phys. Rev. E
65
(
4
),
041906
(
2002
).
9.
Y. L.
Maistrenko
,
B.
Lysyansky
,
C.
Hauptmann
,
O.
Burylko
, and
P. A.
Tass
, “
Multistability in the Kuramoto model with synaptic plasticity
,”
Phys. Rev. E
75
(
6
),
066207
(
2007
).
10.
R.
Berner
,
J.
Fialkowski
,
D.
Kasatkin
,
V.
Nekorkin
,
S.
Yanchuk
, and
E.
Schöll
, “
Hierarchical frequency clusters in adaptive networks of phase oscillators
,”
Chaos
29
,
103134
(
2019
).
11.
R.
Berner
,
J.
Sawicki
, and
E.
Schöll
, “
Birth and stabilization of phase clusters by multiplexing of adaptive networks
,”
Phys. Rev. Lett.
124
(
8
),
088301
(
2020
).
12.
L.
Lücken
,
O. V.
Popovych
,
P. A.
Tass
, and
S.
Yanchuk
, “
Noise-enhanced coupling between two oscillators with long-term plasticity
,”
Phys. Rev. E
93
(
3
),
032210
(
2016
).
13.
T.
Ruangkriengsin
and
M. A.
Porter
, “Low-dimensional analysis of a Kuramoto model with inertia and Hebbian learning,” arXiv:2203.12090 (2022).
14.
S.
Thamizharasan
,
V. K.
Chandrasekar
,
M.
Senthilvelan
,
R.
Berner
,
E.
Schöll
, and
D. V.
Senthilkumar
, “
Exotic states induced by coevolving connection weights and phases in complex networks
,”
Phys. Rev. E
105
(
3
),
034312
(
2022
).
15.
W.
Gerstner
,
R.
Kempter
,
J. L.
van Hemmen
, and
H.
Wagner
, “
A neuronal learning rule for sub-millisecond temporal coding
,”
Nature
383
,
76
(
1996
).
16.
D. C.
Michaels
,
E. P.
Matyas
, and
J.
Jalife
, “
Mechanisms of sinoatrial pacemaker synchronization: A new hypothesis
,”
Circ. Res.
61
(
5
),
704
714
(
1987
).
17.
J.
Buck
and
E.
Buck
, “
Mechanism of rhythmic synchronous flashing of fireflies: Fireflies of southeast asia may use anticipatory time-measuring in synchronizing their flashing
,”
Science
159
(
3821
),
1319
1327
(
1968
).
18.
S. H.
Strogatz
,
D. M.
Abrams
,
A.
McRobie
,
B.
Eckhardt
, and
E.
Ott
, “
Theoretical mechanics: Crowd synchrony on the Millennium Bridge
,”
Nature
438
(
7064
),
43
44
(
2005
).
19.
T.
Danino
,
O.
Mondragón-Palomino
,
L.
Tsimring
, and
J.
Hasty
, “
A synchronized quorum of genetic clocks
,”
Nature
463
(
7279
),
326
330
(
2010
).
20.
C.
Huygens
,
Oeuvres Completes
(
Swets & Zeitlinger Publishers
,
Amsterdam
,
1967
).
21.
E. A.
Martens
,
S.
Thutupalli
,
A.
Fourrière
, and
O.
Hallatschek
, “
Chimera states in mechanical oscillator networks
,”
Proc. Natl. Acad. Sci.
110
(
26
),
10563
10567
(
2013
).
22.
G. H.
Goldsztein
,
L. Q.
English
,
E.
Behta
,
H.
Finder
,
A. N.
Nadeau
, and
S. H.
Strogatz
, “
Coupled metronomes on a moving platform with coulomb friction
,”
Chaos
32
(
4
),
043119
(
2022
).
23.
K.
Wiesenfeld
,
P.
Colet
, and
S.
Strogatz
, “
Frequency locking in Josephson arrays: Connection with the Kuramoto model
,”
Phys. Rev. E
57
(
2
),
1563
1569
(
1998
).
24.
A. F.
Taylor
,
M. R.
Tinsley
,
F.
Wang
,
Z.
Huang
, and
K.
Showalter
, “
Dynamical quorum sensing and synchronization in large populations of chemical oscillators
,”
Science (New York, N.Y.)
323
(
5914
),
614
7
(
2009
).
25.
D.
Călugăru
,
J. F.
Totz
,
E. A.
Martens
, and
H.
Engel
, “
First-order synchronization transition in a large population of relaxation oscillators
,”
Sci. Adv.
6
,
eabb2637
(
2020
).
26.
S.
Danø
,
P. G.
Sørensen
, and
F.
Hynne
, “
Sustained oscillations in living cells
,”
Nature
402
(
6759
),
320
322
(
1999
).
27.
T. M.
Massie
,
B.
Blasius
,
G.
Weithoff
,
U.
Gaedke
, and
G. F.
Fussmann
, “
Cycles, phase synchronization, and entrainment in single-species phytoplankton populations
,”
Proc. Natl. Acad. Sci.
107
(
9
),
4236
4241
(
2010
).
28.
M.
Yamakou
,
P. G.
Hjorth
, and
E. A.
Martens
, “Optimal self-induced stochastic resonance in multiplex neural networks: Electrical versus chemical synapses,”
Front. Comput. Neurosci.
14
,
62
(
2020
).
29.
J.
Klinglmayr
,
C.
Kirst
,
C.
Bettstetter
, and
M.
Timme
, “
Guaranteeing global synchronization in networks with stochastic interactions
,”
New J. Phys.
14
(
7
),
073031
(
2012
).
30.
T.
Gross
and
B.
Blasius
, “
Adaptive Coevolutionary Networks: A Review
,”
J. R. Soc.
5
,
259
271
, (
2008
).
31.
V.
Strauss
, “Stephen Hawking famously said, ‘Intelligence is the ability to adapt to change.’ But did he really say it?” Washington Post, March 29,
2018
.
32.
S. H.
Strogatz
, “
From Kuramoto to Crawford: Exploring the onset of synchronization in populations of coupled oscillators
,”
Physica D
143
,
1
20
(
2000
).
33.
M. A.
Gkogkas
,
B.
Jüttner
,
C.
Kuehn
, and
E. A.
Martens
, “
Graphop mean-field limits and synchronization for the stochastic Kuramoto model
,”
Chaos
32
,
113120
(
2022
).
34.
E. A.
Martens
and
K.
Klemm
, “
Transitions from trees to cycles in adaptive flow networks
,”
Front. Phys.
5
,
1
10
(
2017
).
35.
E. A.
Martens
and
K.
Klemm
, “Cyclic structure induced by load fluctuations in adaptive transportation networks,” in Progress in Industrial Mathematics at ECMI 2018, edited by I. Faragó, F. Izsák, and P. L. Simon (Springer International Publishing, 2019), pp. 147–155.
36.
H.
Mestre
,
T.
Du
,
A. M.
Sweeney
,
G.
Liu
,
A. J.
Samson
,
W.
Peng
,
K. N.
Mortensen
,
F. F.
Stæger
,
P. A. R.
Bork
,
L.
Bashford
,
E. R.
Toro
,
J.
Tithof
,
D. H.
Kelley
,
J. H.
Thomas
,
P. G.
Hjorth
,
E. A.
Martens
, and
R. I.
Mehta
, “
Cerebrospinal fluid influx drives acute ischemic tissue swelling
,”
Science (New York, N.Y.)
367
(
6483
),
eaax7171
(
2020
).
37.
T.
Bohr
,
P. G.
Hjorth
,
S. C.
Holst
,
S.
Hrabětová
,
V.
Kiviniemi
,
T.
Lilius
,
I.
Lundgaard
,
K.-A.
Mardal
,
E. A.
Martens
,
Y.
Mori
,
U. V.
Nägerl
,
C.
Nicholson
,
A.
Tannenbaum
,
J. H.
Thomas
,
J.
Tithof
,
H.
Benveniste
,
J. J.
Iliff
,
D. G.
Kelley
, and
M.
Nedergaard
, “
The glymphatic system: Current understanding and modeling
,”
iScience
25
(
9
),
104987
(
2022
).
38.
J. P.
Taylor-King
,
D.
Basanta
,
S. J.
Chapman
, and
M. A.
Porter
, “
Mean-field approach to evolving spatial networks, with an application to osteocyte network formation
,”
Phys. Rev. E
96
(
1
),
012301
(
2017
).
39.
B.
Skyrms
and
R.
Pemantle
, “A dynamic model of social network formation,” in Adaptive networks (Springer, 2009), pp. 231–251.
40.
W.
Gerstner
and
W. M.
Kistler
, “
Mathematical formulations of Hebbian learning
,”
Biol. Cybern.
87
(
5–6
),
404
415
(
2002
).
41.
A.
Goriely
,
E.
Kuhl
, and
C.
Bick
, “
Neuronal oscillations on evolving networks: Dynamics, damage, degradation, decline, dementia, and death
,”
Phys. Rev. Lett.
125
(
12
),
128102
(
2020
).
42.
N.
Caporale
and
Y.
Dan
, “
Spike timing-dependent plasticity: A Hebbian learning rule
,”
Annu. Rev. Neurosci.
31
,
25
46
(
2008
).
43.
B.
Duchet
,
C.
Bick
, and
A.
Byrne
, “Mean-field approximations with adaptive coupling for networks with spike-timing-dependent plasticity,” bioRxiv (2022), p. 2022-07.
44.
Note that this property is due to the model’s special form that is invariant with regard to stretching time; other models, e.g., with intrinsic dynamics of the form ϕ ˙ l = f ( ϕ l ) + l g ( ϕ l , ϕ m ) do not necessarily share this property, an example being the Theta neuron model.
45.
Y.
Kuramoto
,
Chemical Oscillations, Waves, and Turbulence
(
Springer-Verlag
,
New York
,
1984
).
46.
A.
Dhooge
,
W.
Govaerts
,
Y. A.
Kuznetsov
,
H. G. E.
Meijer
, and
B.
Sautois
, “
New features of the software MatCont for bifurcation analysis of dynamical systems
,”
Math. Comput. Modell. Dyn. Syst.
14
,
147
175
(
2008
).
47.
C.
Kuehn
,
Multiple Time Scale Dynamics
(
Springer
,
2015
), Vol. 191.
48.
P.
Ashwin
and
O.
Burylko
, “
Weak chimeras in minimal networks of coupled phase oscillators
,”
Chaos
25
(
1
),
013106
(
2015
).
49.
B.
Jüttner
,
C.
Henriksen
, and
E. A.
Martens
, “
Birth and destruction of collective oscillations in a network of two populations of coupled type 1 neurons
,”
Chaos
31
,
023141
(
2021
).
50.
M.
Thiele
,
R.
Berner
,
P. A.
Tass
,
E.
Schöll
, and
S.
Yanchuk
, “Asymmetric adaptivity induces recurrent synchronization in complex networks,”
Chaos
33
,
023123
(
2023
).
51.
D.
Kasatkin
and
V.
Nekorkin
, “
Dynamics of the phase oscillators with plastic couplings
,”
Radiophys. Quantum Electron.
58
(
11
),
877
891
(
2016
).
52.
O.
Burylko
,
E. A.
Martens
, and
C.
Bick
, “
Symmetry breaking yields chimeras in two small populations of Kuramoto-type oscillators
,”
Chaos
32
(
9
),
093109
(
2022
).
53.
D. M.
Abrams
,
R.
Mirollo
,
S. H.
Strogatz
, and
D. A.
Wiley
, “
Solvable model for chimera states of coupled oscillators
,”
Phys. Rev. Lett.
101
,
084103
(
2008
).
54.
M. J.
Panaggio
,
D. M.
Abrams
,
P.
Ashwin
, and
C. R.
Laing
, “
Chimera states in networks of phase oscillators: The case of two small populations
,”
Phys. Rev. E
93
(
1
),
012218
(
2016
).
55.
E. A.
Martens
,
C.
Bick
, and
M. J.
Panaggio
, “
Chimera states in two populations with heterogeneous phase-lag
,”
Chaos
26
(
9
),
094819
(
2016
).
56.
H.
Hong
and
S. H.
Strogatz
, “
Kuramoto model of coupled oscillators with positive and negative coupling parameters: An example of conformist and contrarian oscillators
,”
Phys. Rev. Lett.
106
(
5
),
054102
(
2011
).
57.
C.
Bick
,
M. J.
Panaggio
, and
E. A.
Martens
, “
Chaos in Kuramoto oscillator networks
,”
Chaos
28
,
071102
(
2018
).
58.
J. D.
Crawford
, “
Amplitude expansions for instabilities in populations of globally-coupled oscillators
,”
J. Stat. Phys.
74
(
5
),
1047
1084
(
1994
).
59.
E.
Andreas Martens
,
E.
Barreto
,
S. H.
Strogatz
,
E.
Ott
,
P.
So
, and
T. M.
Antonsen
, “
Exact results for the Kuramoto model with a bimodal frequency distribution
,”
Phys. Rev. E
79
(
2
),
026204
(
2009
).
60.
D.
Pazó
and
E.
Montbrió
, “
Existence of hysteresis in the Kuramoto model with bimodal frequency distributions
,”
Phys. Rev. E
80
(
4
),
046215
(
2009
).
61.
B.
Pietras
,
N.
Deschle
, and
A.
Daffertshofer
, “
Equivalence of coupled networks and networks with multimodal frequency distributions: Conditions for the bimodal and trimodal case
,”
Phys. Rev. E
94
(
5
),
052211
(
2016
).
62.
H.
Hong
and
E. A.
Martens
, “
A two-frequency-two-coupling model of coupled oscillators
,”
Chaos
31
(
8
),
083124
(
2021
).
63.
M. A.
Gkogkas
,
C.
Kuehn
, and
C.
Xu
, “Mean field limits of co-evolutionary heterogeneous networks,” arXiv:2202.01742 (2022).
64.
M. A.
Gkogkas
,
C.
Kuehn
, and
C.
Xu
, “Continuum limits for adaptive network dynamics,” arXiv:2109.05898 (2021).
65.
C.
Kuehn
and
C.
Xu
, “Vlasov equations on digraph measures,” arXiv:2107.08419 (2021).
66.
J.
Fialkowski
,
S.
Yanchuk
,
I. M.
Sokolov
,
E.
Schöll
,
G. A.
Gottwald
, and
R.
Berner
, “
Heterogeneous nucleation in finite-size adaptive dynamical networks
,”
Phys. Rev. Lett.
130
(
6
),
067402
(
2023
).
67.
C.
Bick
and
D.
Sclosa
, “Mean-field limits of phase oscillator networks and their symmetries,” arXiv:2110.13686 (2021).