Asymmetric Adaptivity induces Recurrent Synchronization in Complex Networks

Rhythmic activities that alternate between coherent and incoherent phases are ubiquitous in chemical, ecological, climate, or neural systems. Despite their importance, general mechanisms for their emergence are little understood. In order to ﬁll this gap, we present a framework for describing the emergence of recurrent synchronization in complex networks with adaptive interactions. This phenomenon is manifested at the macroscopic level by temporal episodes of coherent and incoherent dynamics that alternate recurrently. At the same time, the dynamics of the individual nodes do not change qualitatively. We identify asymmetric adaptation rules and temporal separation between the adaptation and the dynamics of individual nodes as key features for the emergence of recurrent synchronization. Our re-sults suggest that asymmetric adaptation might be a fundamental ingredient for recurrent synchronization phenomena as seen in pattern generators, e.g., in neuronal systems.


I. INTRODUCTION
2][3][4][5][6] Structural changes in such networks depend on the state of individual nodes.The state of the nodes in turn depends on the connectivity.This results in a feedback loop between the dynamics (function) and the network structure.In neuronal networks, for example, spike-timing-dependent plasticity, a special type of network adaptation, [7][8][9][10][11][12][13] contributes to learning 14 or antikindling 15 effects.Another adaptation mechanism, structural plasticity, is responsible for a homeostatic regulation of electrical activity in the brain. 16Beyond neural networks, adaptive networks are used in communication systems 17 and for modeling complex behavior in social, climate, or ecological systems. 5,18,19n realistic systems, the adaptation rules are likely to be heterogeneous.For example, the adaptation rule in neuronal networks may depend on the distance between interacting neurons. 20,21So far, little is known about the dynamical effects induced by heterogeneous adaptation. 22,23The question arises as to which extent heterogeneous adaptation produces new functionality and, thus, represents an important element for system behavior.In this work, we answer this question by showing that heterogeneous adaptation can be a determining ingredient for producing new collective network function.In particular, we describe how heterogeneous adaptation induces the phenomenon of recurrent synchronization.
Recurrent synchronization, which we study here, is a macroscopic phenomenon in dynamical networks involving the recurrent switching between synchronous behavior, represented, e.g., by phase-locking, and asynchronous behavior, represented, e.g., by frequency clustering.During recurrent synchronization, a macroscopic observable exhibits bursting behavior, whereas the individual nodes are not required to burst at the microscopic level.In our study, we consider the nodes to have simple oscillatory dynamics.Therefore, recurrent synchronization as a macroscopic effect contrasts bursting found in neuronal networks, [24][25][26][27] where single neurons can have alternating periods of quiescence and fast spiking.
Other studies [28][29][30][31][32] have reported the occurrence of a similar phenomenon, called collective bursting, in neuronal networks.An important difference of Refs.28-32 from our work, however, is that in these studies collective bursting phenomenon is induced and observable on the microscopic level of individual neurons, whereas recurrent synchronization is not manifested by bursting on the microscopic level.Another adaptation-induced switching between phase-locking and periodic oscillation has been observed in the work of Ref. 33, which is related to fold-homoclinic bursting. 34In a two-population excitatory-inhibitory network of quadratic integrate and fire neurons, a dynamical phenomenon was observed showing bursts of high frequency and high amplitude activity at a slow burst rate. 35In contrast to recurrent synchronization considered here, the populations change their individual dynamics from oscillatory to steady state during the episodes of high activity and quiescence.This phenomenon of activity bursting has been also reported in a population of excitable units adaptively coupled to a pool of resources. 36n a broader sense, the mechanism for recurrent synchronization, which we report here, is based on a recurrent slow dynamics of hidden variables that we relate to adaptive coupling weights between dynamical populations in this work.When adaptation is heterogeneous (asymmetric), the hidden variables lead to reemergence of episodes of synchronization and desynchronization.Effects such as converse symmetry breaking 37,38 or asymmetryinduced synchronization 39 have been described only very recently and are about to change our understanding of "good synchronization conditions" for dynamical systems.In this light, recurrent synchronization adds a new layer of dynamical complexity that can be induced by asymmetry and, hence, heterogeneity, in complex dynamical networks.
In our study, we use models with different hierarchical organization.While two neurons with adaptive coupling is a sufficient minimal model, other models show that recurrent synchronization can occur in dynamic models of varying complexity.A rather complex system consists of two populations of Hodgkin-Huxley neurons with different adaptation rules within and between populations.
The adaptation rules are given by spike-timing-dependent plasticity with different adaptation functions.We also propose reduced phenomenological models of two coupled Hodgkin-Huxley neurons as well as phase oscillators equipped with slowly evolving coupling weights.We show that the reduced models are capable of describing the main dynamical features arising in the intra-and inter-population dynamics.While the more complex model is studied numerically, the reduced phenomenological models are analyzed in more detail by methods from geometric singular perturbation, averaging, and bifurcation theories.In this work, we propose a general methodology to study recurrent synchronization in systems with arbitrary complex dynamical nodes and continuous as well as noncontinuous adaptation rules.
Our study reveals collective mechanisms behind the appearance of recurrent synchronization.In particular, we explain the importance of the following ingredients for the emergence of recurrent synchronization: slow adaptation, i.e., the timescale separation between the adaptation and the individual neuronal dynamics; asymmetry of adaptation rules; recurrent (periodic/spiking) dynamics of individual neurons.We hypothesize that asymmetric adaptivity might play a fundamental role in emergence and impairment of neuronal pattern generators.

A. Network of Hodgkin-Huxley neurons
The dynamics for the network of synaptically coupled Hodgkin-Huxley neurons are given by the following equations: 40,41 Here, V i is the membrane potential of the ith neuron, and E Na = 50 mV, E K = −77 mV, and E L = −54.4mV are reverse potentials of different ion channels.The dynamic variables m i , h i , and n i denote gating variables describing the probability of open ion-channels for sodium activation, sodium deactivation, and potassium, respectively.The membrane capacity is given by C = 1 µF/cm 2 .We set the reverse potential of coupling E r = 20 mV, which corresponds to excitatory neurons.The functions α i and β i are given as 65  20   , .
The corresponding conductance parameters are set to g Na = 120 mS/cm 2 , g K = 36 mS/cm 2 , and g L = 0.3 mS/cm 2 .The constant input currents are randomly chosen from the intervals [4.99 µA, 5.01 µA] and [12.99 µA, 13.01 µA] for populations 1 and 2, respectively.Different input currents I i account for different spiking rates of the two populations and different neurons.
The coupling terms include the synaptic activity s j of the presynaptic neuron, multiplied by the coupling weights κ ij , which changes discontinuously due to synaptic plasticity every time the ith or jth neuron exhibits a spike.Whenever a neuron spikes, the coupling weights between the neurons receive an update according to κ ij → κ ij + δ • W ( T) (δ = 0.005), which depends on the timing difference T ij = t i − t j between the last spiking events t i and t j of the ith and jth neuron, respectively.We identify a spiking event by checking whether the voltage of the neuron in question exceeds 0 in the given time step.For the two populations, we use different update functions W. For population 1 (coupling weights of all connections with postsynaptic neurons in population 1), the update function is anti-Hebbian 42 and of the following form: ( The parameters used for Figs. 2 and 6 are A 1 = 1.17,A 2 = 0.4, τ 1 = 0.25, τ 2 = 1.1, and e = exp(1).For population 2, the update function is symmetric 43 and of the following form: with c p = 1.5, c d = 0.53, τ p = 1.8, and τ d = 5.To ensure bounded coupling weights, we restrict them to the interval [0, 1.5].The two plasticity rules exemplify the diversity of adaptation rules known for nervous systems.We note, however, that the dynamical system introduced in this work can be conceived as a toy model possessing complex dynamical features rather than a model of a realistic system.

B. Two phase oscillators
As each population in the system described above is likely to synchronize, we also consider a simpler paradigmatic model for the dynamics of the synchronous population: the phase oscillator. 44,45aking into account the adaptation, the model of two adaptively coupled phase oscillators reads 46,47 where φ i ∈ [0, 2π) represents the phases of the oscillatory dynamics for each population.When the phase increases from 0 to 2π , this corresponds to one oscillation period.The natural frequencies of the ith oscillator are denoted by ω i .The adaptive coupling weights are given by κ i .The phase lag parameter α is included that may account for a small information transmission delay. 48The parameters a and b scale the influence of the phase dynamics on the coupling weights.An additional asymmetry in the adaption rules for κ 1 and κ 2 is introduced by the parameter β.

C. Macroscopic observables
To measure the collective motion of the whole network, we introduce an observable that provides an average of the individual node dynamics.The behavior of each dynamical node x i (t) can be mapped to a motion on the unit circle represented by a phase variable φ i (t). 44Particularly, for Hodgkin-Huxley neurons considered in Sec.II A, each node possesses a sequence of spiking events {t i,1 , t i,2 , . . .}.The mapping used in this study is given by . The collective observable is then defined by the Kuramoto order parameter If, for a certain time t, the phases are incoherently spread over the interval [0, 2π), the order parameter is zero while it is unity if all phases are the same.Furthermore, by looking at the temporal behavior of R(t), we are able to distinguish between phase-locking, which corresponds to an order parameter approximately constant in time, i.e., R(t) ≈ const, and the loss of phase relation that corresponds to an order parameter oscillating between incoherence (small values of R) and coherence R ≈ 1.
For the network of Hodgkin-Huxley neurons, we consider the firing density as a second macroscopic observable.The firing density is determined by dividing the time period of observation into time bins of width = 3000 ms and counting the spikes for each time bin.Afterward, the count per time bin is divided by the number of neurons in the group.

III. RECURRENT SYNCHRONIZATION
In this section, we introduce the phenomenon of recurrent synchronization, where certain macroscopic observables display bursting behavior while the individual nodes emit single periodic spikes.We show that such a behavior is caused by alternating episodes of synchronization and desynchronization in a dynamical network.
The reported phenomenon can be described in general terms, omitting nonessential details of specific model implementation.The corresponding Fig. 1 summarizes the main phenomenology.We distinguish between the microscopic (individual) and the macroscopic (collective) scales.While the individual nodes show periodic oscillatory behavior [Fig.1(b)], the collective motion exhibits an alternation between time intervals of high activity, i.e., fast changes of the collective observable, and time intervals of low activity [Fig.1(c)], i.e., comparatively slow changes of the collective observable.
The mechanism for the recurrent intervals of synchronization and desynchronization can be explained by additional, hidden, slow variables [Fig.1(d)].Such hidden variables are not necessarily observable, but they play the key role as the internal control variables that determine the synchronization level.As a result of a subtle interplay between the nodal dynamics and the hidden variables, alternating transitions between synchronized (phase-locked state) and desynchronized episodes (bursting state) can be induced [Fig.1(d)].As described below, the slow hidden variables naturally arise in networks with slow network adaptivity.
The recurrent synchronization phenomenon, which we report here, differs from those observed in Refs.28-32, 34, 49, 50, where single neurons possess alternating periods of spiking and rest.In our case, the nodes exhibit periodic behavior (tonic spiking) at all times, which makes the observed phenomenon even more surprising.

A. Recurrent synchronization in interacting populations of Hodgkin-Huxley neurons with spike-timing-dependent plasticity
Models of interacting populations are well-known paradigms for studying dynamics in many real-world systems on the mesoscopic scale. 51,52These models have particular importance in neuroscience for the modeling of interacting brain regions or other functional units 53 and are also used in social sciences with autonomous agents interacting with each other based on their population affiliation. 54,55o show recurrent synchronization in a complex dynamical network setup, we implement two populations of Hodgkin-Huxley neurons with spike-timing-dependent plasticity (STDP) (Fig. 2), where the two populations differ in the adaptation functions of the coupling weights and in the input current, as described in Sec.II A and shown in Fig. 2(a).
Figures 2(c) and 2(e) show the temporal behavior of two macroscopic variables: the order parameter R(t) that measures the level of synchronization and the firing density measuring the mean level of activity of the whole system and the individual populations, respectively.The occurrence of several transitions between a nearly constant (phase-locking episode) to a strongly oscillating order parameter (bursting episode) is clearly visible.Figure 2(d) presents a zoom into a bursting episode, showing many oscillations in a small time interval of 1 s. Figure 2(e) also displays the firing densities for two different populations showing strong differences for the episodes of a fast oscillating order parameter.Figures 2(f) and 2(g) depict raster plots for the two dynamical episodes.It is clear that the collective bursting phenomenon in the order parameter is not due to the bursting behavior of individual neurons that show consistent tonic spiking at any time.The phenomenon is caused by frequency desynchronization between the two populations.The latter can be seen in the different inter-spike intervals of the populations.These results are also robust to noisy inputs, which we model with random α-spikes, and higher heterogeneity in the input currents as well as heterogeneity in the update functions of the coupling weights, see Sec.III B for details.
Interestingly, a sharp drop in spiking activity coincides with both the onset and the termination of the bursting phase.

B. Robustness of recurrent synchronization in networks of Hodgkin-Huxley neurons
To investigate the robustness of recurrent synchronization in two populations of Hodgkin-Huxley neurons, we first simulate a network of N = 200 neurons where every neuron receives an independent input in the voltage dynamics.We choose an α-train input of the form 56 modeling random post-synaptic potentials arriving at the ith neuron with the time intervals between two potentials being distributed as N (14 ms, 4 ms).The time evolution of the order parameter for I = 0.01 can be seen in Fig. 3. Recurrent synchronization is still present, but the oscillations in the phase-locked periods have larger amplitudes than in the noiseless case shown in Fig. 2(b).Furthermore, the results shown in Fig. 2(c) are also present for increased heterogeneity in the input current each neuron receives, which is shown in Fig. 4. Here, we show the order parameter of the network for the last 300 s of a 600 s long simulation for non-homogeneous input currents uniformly drawn from the intervals [4.5 µA, 5.5 µA] and [12.5 µA, 13.5 µA], for population 1 and 2, respectively.Recurrent synchronization is still present in this more heterogeneous system and for the longer time interval.Recurrent synchronization is likewise existent for heterogeneity in the update functions of the coupling weights.In Fig. 5, the order parameter for a network of 200 neurons with distributed update functions is shown.Recurrent synchronization is again observable even for the case of distributed input currents, update functions, and initial conditions, demonstrating the overall robustness of recurrent synchronization in heterogeneous setups.

C. Phenomenological mean field reduction
The simulation in Fig. 2 indicates that the dynamics within the populations is rather coherent while it can exhibit incoherence between the populations.In order to deepen the understanding of the emergence of recurrent synchronization, we consider the reduced model of the following general form: where the mean field variables x 1 and x 2 describe the coherent dynamics within the populations 1 and 2, respectively.The variable coupling weights κ 1 and κ 2 are the coupling weights representing the mean inter-population connections.We consider a small parameter 0 < 1, which enables the separation of time scales between the fast dynamics of the nodes and the slow dynamics of the coupling weights.The adaptation dynamics is assumed to be on the same slow timescale 1/ .The functions F 1,2 govern the dynamical behavior of the corresponding populations and W 1,2 determine the adaptation rules of the coupling weights.
Using the explicit splitting of the time scales (smallness of ), we develop a further reduction in the mean-field equations ( 10)- (13).In particular, as the coupling weights (κ 1 , κ 2 ) evolve slowly, we consider them as parameters for the dynamics of the oscillators.Depending on the coupling weights, the coupled oscillator system ẋ1 = F 1 (x 1 , x 2 , κ 1 ), ẋ2 = F 2 (x 2 , x 1 , κ 2 ) shows either synchronized or desynchronized motion.
Denoting by • the temporal averaging on the time scale that is much longer compared to the fast changes of the variables x 1,2 (t) but shorter than the timescale 1/ of the adaptation, we obtain and κ1,2 ≈ κ1,2 .This leads to the following two-dimensional reduced model: describing the dynamics of the coupling weights κ 1 and κ 2 on the slow scale.In system ( 14), we have additionally rescaled time t s = t and, thus, removed the factors .The prime in κ 1 and κ 2 denotes the derivative with respect to the new slow time t s .The coupling weights κ 1 and κ 2 play the role of hidden variables governing the transitions between synchronization and desynchronization as shown schematically in Fig. 1(d).System ( 14) allows for a detailed bifurcation analysis of the mechanism behind the emergence of recurrent synchronization.Clearly, for recurrent synchronization, the dynamics of the hidden variables must exhibit recurrent motions.As shown in the following, such a recurrent motion can be observed when the adaptation functions between the two dynamical populations are non-symmetric.
As the system (14) describes the (hidden) dynamics with respect to the slow time, it does not involve the small parameter and, thus, does not involve this time scale separation and can be efficiently computed.
In Secs.IV and V, we derive the models for slowly evolving coupling weights for two specific examples: two adaptively coupled phase oscillators and the more complex case of two coupled Hodgkin-Huxley neurons equipped with spike timing-dependent plasticity.The latter application shows, moreover, the generality of our approach with regards to complex local dynamics and adaptation rules.
From the mathematical point of view, the reduction is based on geometric singular perturbation theory 57,58 and averaging 59 for synchronized and desynchronized regimes, respectively.Both theoretical approaches are implemented analytically as well as semianalytically.
The averaging of the local dynamics is easily made for the case when the dynamics of the subsystems ( 10) and ( 11) (with fixed κ 1 and κ 2 ) converge to an equilibrium (x * 1 (κ 1 , κ 2 ), x * 2 (κ 1 , κ 2 )).The procedure is known as "adiabatic elimination," and subsystems ( 10) and ( 11) as the "layer equation" or "fast subsystem." 57If the fast subsystem has one or more equilibria, then they satisfy the following equations: where the union of all solutions (x * 1 , x * 2 ) compose the so-called critical manifold.If an equilibrium (x * 1 (κ 1 , κ 2 ), x * 2 (κ 1 , κ 2 )) is stable with respect to the fast subsystem (10) and (11), then the solutions will rapidly converge to this equilibrium, and the reduced system (14)   can be simply obtained by substituting this equilibrium point into (12) and ( 13), If the fast subsystem ( 10) and (11) does not converge to a stable equilibrium, the adiabatic elimination cannot be applied.In such a case, direct averaging must be used.If, for fixed (κ 1 , κ 2 ), the solution to (10) and ( 11) converges to a periodic state x 1 (t; κ 1 , κ 2 ), x 2 (t; κ 1 , κ 2 ) with a period T(κ 1 , κ 2 ), then this averaging procedure leads to

IV. REDUCED MODEL: TWO HODGKIN-HUXLEY NEURONS WITH SPIKE-TIMING-DEPENDENT PLASTICITY
When the dynamical populations in Fig. 2 are synchronized, the collective dynamics of each population can be approximated by a single Hodgkin-Huxley neuron.As a result, the dynamical network can be reduced to two coupled Hodgkin-Huxley neurons with asymmetric synaptic plasticity.The synaptic weight defined as κ 1 = κ 12 is updated accordingly to the STDP rule with the update function W 1 and κ 2 = κ 21 with W 2 (see Sec. II A), as in the case of two populations in Fig. 2.
For such a system of two coupled Hodgkin-Huxley neurons with STDP, the reduction procedure from Sec. III C can be applied numerically by using the averaging method in (17) for sufficiently large T.
Practically, the slow flow for the two Hodgkin-Huxley neurons in Fig. 6(d) is obtained using a grid in the plane of the coupling weights (κ 1 , κ 2 ) and estimating the averages in (17) for each grid point.For this, the distribution ρ(ξ ) of inter-spike intervals T is first computed.Then, the right-hand side of the reduced systems are estimated as Note that any scaling factors in (18), including δ, can be ignored as they can be scaled out by time rescaling, similar to the rescaling of the small parameter .Importantly, the fast systems have to be computed once in order to find the distributions ρ(ξ ) for every grid point in the κ 1 , κ 2 plane.Once these data are created, the slow flow can be computed using (18) for any update function without simulating the fast system.This makes the study of the effects of different plasticity rules much more feasible.
Figure 6(a) displays recurrent synchronization with a zoom into an asynchronous episode shown in Fig. 6(b).The reduced flow for the coupling weights ( 14) is presented in Fig. 6(d) along with the projection of the simulated trajectory onto the (κ 1 , κ 2 )-plane.The recurrent synchronization, found in simulations, corresponds to a periodic trajectory in the reduced model (green line).The trajectory in Fig. 6(d) enters synchronous and asynchronous regimes recurrently.Moreover, the averaged flow (black lines) and the projected flow based on the simulation of the whole system are in very good agreement, and only small deviations near the boundary can be observed.Therefore, our proposed reduction method that eliminates the fast dynamics of the oscillators while keeping the slow dynamics of the coupling weights can also be applied to discontinuous adaptation functions.
A closer view into transition from synchronization to desynchronization can be seen in Fig. 6(c).At the end of a synchronized episode, the oscillations of the order parameter are gradually building up from small to large variations.This observation could be a useful indicator for the onset of bursting events where small variations of the order parameter may be interpreted as precursors.
A numerical bifurcation diagram for parameters τ 1 (parameter of W 1 ) and γ (parameter of W 2 ) is presented in Fig. 6 The regions G1 and G2 (blue) mark areas where the equilibrium (κ 1 , κ 2 ) = (0, 0) is stable; C (red): a stable equilibrium on the critical manifold; A (orange): a stable recurrent synchronization limit cycle; E (green): a stable limit cycle on the critical manifold.The other regions mark areas of bistability of the equilibria described above.Also shown are analytically obtained bifurcation lines for the whole system (red) and the (0,0) equilibrium in the reduced system (blue), with SN denoting a Saddle-Node bifurcation, PF a Pitchfork bifurcation and H a Hopf bifurcation (dashed lines).The boundaries for regions A, B, E, and F are obtained numerically.The diagrams in the two columns on the right and the row in the bottom show exemplary trajectories for the reduced system in the (κ 1 , κ 2 ) plane for the regions denoted in the bifurcation diagram.Black lines mark stable trajectories, blue lines mark unstable trajectories, and green lines mark the limit cycle.The (0,0) equilibrium can either be stable (red and filled) or unstable (red and void).The blue areas denote the basins of attraction of the (0,0) equilibrium in the case of bistability.Green dots indicate stable fixed points on the critical manifold and the red line marks the boundary of the critical manifold.
recurrent synchronization is observed for a wide range of parameters.The periodic motion in the hidden variables shown in Fig. 6(d) is induced by asymmetric adaptation [see Fig. 2(b)] and is not be found for symmetric adaptation rules.For comparison, Fig. 7 illustrates how the system of two Hodgkin-Huxley neurons evolves when the update functions for κ 1 and κ 2 are identical.For the considered parameter values, we observe no recurrent synchronization.Shown is the trajectory of the whole system (in green) projected onto the (κ 1 , κ 2 )-plane with both update functions being equal to W 1 defined in Eq. ( 2) [Fig.7(a)] and to W 2 defined in Eq. ( 3) [Fig.7(b)].These trajectories are in very good agreement with the reduced flow based on the averaging approach described at the beginning of this section (black lines with arrows).Both the trajectory and the flow show no presence of a periodic orbit traversing the boundary between synchronous and asynchronous spiking regimes, which would be the way recurrent synchronization appears in this representation.
Figure 7 also displays exemplary distributions of the spike timing differences: Figs.7(c) and 7(d) for asynchronous regimes, Fig. 7(e) near the boundary, and Fig. 7(f) for the synchronous regime, with the two peaks consistent with phase-locked spiking.These figures illustrate the impact the coupling weights have on the relative spiking times.

V. REDUCED MODEL: TWO PHASE OSCILLATORS
Applying the averaging procedure from Sec. III C to two coupled phase oscillators (4)- (7), we obtain the following twodimensional system: where are the corresponding temporal averages of adaptation functions.We emphasize that the explicit dependence of the functions W 1,2 on the coupling weights stem from the dependence of the fast dynamics φ 1 (t) and φ 2 (t) on these weights.The functions W 1,2 (κ 1 , κ 2 ) are computed explicitly in Appendices A and B, which allow the properties of the reduced system (19) to be studied in detail.While the technical details of the reduction are described in Secs.A and B, we note here that this procedure is performed by two different methods: adiabatic elimination for such parameter values κ 1 , κ 2 that the subsystem ( 4) and ( 5) is synchronized (phase-locked) and the averaging along a periodic orbit when ( 4) and ( 5) is not synchronized.
The bifurcation diagram for system (19) in Fig. 8 shows the regions for which recurrent synchronization is observed (yellow) in the (a, b)-parameter plane for fixed parameters ω = ω 1 − ω 2 and α.In the parameter region A, the recurrent synchronization is the only stable regime [Fig.8(b)], while in region B, recurrent synchronization coexists with an effectively uncoupled state, i.e., equilibrium at (κ 1 , κ 2 ) = (0, 0) of the reduced system [Fig.8(d)].
Figures 8(c) and 8(e) show the time evolution of the order parameter (8) in the regime of recurrent synchronization.Interestingly, we observe the repeated occurrence of two synchronous regimes: an in-phase and an anti-phase locking interrupted with the running phase regime.In Fig. 8(f), yet another stable state (red line) is shown corresponding to desynchronization of the two oscillators leading to an effective decoupling.The corresponding states in the reduced model are displayed in the phase portraits in Figs.8(b) and 8(d).These phase portraits of (19) show limit cycles passing through regions of synchronous and asynchronous relative phase dynamics, which explains the dynamics of the order parameter in Figs.8(c) and 8(e).
Recurrent synchronization can also be found for other values of the asymmetry parameter β, see Fig. 12, and time-scale separation parameter , see Fig. 11.In particular, the symmetric case β = 0 reveals no recurrent synchronization.It seems that a certain amount of asymmetry in the adaptation functions is necessary to obtain a stable recurrent synchronization as the regions for recurrent synchronization lie off the lines |a| = |b| in Fig. 8(a).
The bifurcation analysis reveals that the emergence of recurrent synchronization is the result of a fold bifurcation of the limit

ARTICLE
scitation.org/journal/chacycle of the reduced system, as can be seen in Fig. 9, which shows a much larger area in the (a, b)-parameter plane.Here, all regions with distinct stable states are displayed with exemplary trajectories for the reduced system.Furthermore, analytical bifurcation lines are drawn in red and blue, with the dashed lines corresponding to Hopf bifurcations and the solid lines to Saddle-Node-or Pitchfork bifurcations.The derivation of these lines can be found in Appendix C. The diagram describing the border between regions G2 and B shows the trajectories after the annihilation of the two limit cycles with only the (κ 1 , κ 2 ) = (0, 0) point remaining as a stable state.Recurrent synchronization can also appear, when the synchronous relative phase dynamics become unstable.This can either happen, when a "slow" limit cycle, i.e., a limit cycle with changes in the relative phase dynamics on the order of , crosses the boundary between the synchronous and asynchronous regime (transitions from region E to A and F to B, respectively) or via the destabilization of a stable fixed point at the boundary between the two regimes.
To check whether the recurrent synchronization limit cycle contains additional information about the bifurcation present in the system, the periods are shown in Fig. 10.Here, the periods are colorcoded in the parameter plane (a.b), with the light areas signifying no recurrent synchronization.We can identify higher periods for lower values of a and b and vice versa.Since a and b denote the amplitude of the phase influence on the dynamics of the coupling weights, they, therefore, impact the speed with which trajectories move in the reduced phase space of the coupling weights (κ 1 , κ 2 ).This yields an inverse relationship between the values of a, b and the period of the limit cycle.

VI. CONCLUSION
In this work, we have described recurrent synchronization for two populations of adaptively coupled oscillators.The phenomenon is shown to be robust to noise, variation of parameters, initial states, and heterogeneity of oscillators.
To shed light on the emergence of recurrent synchronization, we have introduced a mean-field model approach.A detailed bifurcation analysis of a simplified phenomenological phase oscillator model provides remarkable insights into the mechanisms leading to recurrent synchronization.The separation between the time scales of the adaptation and the fast individual dynamics allows for the model reduction in the case of two phase oscillators and two Hodgkin-Huxley neurons.The proposed reduction approach is applicable to systems beyond phase oscillators and, in this regard, it methodologically generalizes recently used methods from Ratas et al. 60 and Franović et al. 61 Our study is further complemented by an analysis of a mean-field model consisting of coupled Hodgkin-Huxley neurons equipped with spike-timing-dependent plasticity.For this more complicated system, our numerical averaging approach makes the study of the effects of plasticity on neuronal dynamics much more feasible.
Heterogeneity in the adaptation rule has been known to exist, e.g., in neural systems for a long time, 62 but their dynamical effects are only rarely studied and analytical methods are widely unexplored. 22,23ur findings contribute to the important question how heterogeneity may change the dynamics of complex networks.We have found that such heterogeneity in the adaptivity is able to significantly change the macroscopic dynamics of a dynamical network.Moreover, our analysis contributes to the identification of the onset of critical phenomena and extreme events.In the transition from synchrony to asynchrony, we have observed gradually growing oscillations of the order parameter.This observation may be used as characterizing feature for the onset of the bursting period and, thus, may indicate substantial changes in the properties of the network.
Our results may contribute to the understanding of pattern generators, their disease-related impairment, and, specifically, the origin of Parkinsonian resting tremor.4][65] Neuronal discharges in the tremor frequency range in different parts of the basal ganglia as well as the thalamus were found to be coherent and transiently locked to the peripheral tremor. 66Furthermore, causal data analysis suggests that this activity drives the muscular tremor, as opposed to being just driven by the sensory feedback from the periphery. 67In functional magnetic resonance imaging (fMRI) studies, it was shown that basal ganglia get active transiently when tremor episodes emerge, whereas activity in the cerebello-thalamo-cortical circuit is tightly related to the tremor amplitude. 65According to the "dimmer-switch" hypothesis, the basal ganglia activate the tremor (like a light switch) and the cerebello-thalamo-cortical circuit modulates tremor amplitude (like a light dimmer). 65,68However, it remains elusive what causes spontaneous switching between epochs with alternating (anti-phase) tremor and synchronous (in-phase) tremor, connected by epochs without stable phase relationship. 63,69Based on our results, we suggest that asymmetric adaptivity involving two populations (Fig. 2) responsible for two antagonistic muscles might cause the clinically observed recurrent epochs of synchrony with different phase differences, e.g., in-phase and anti-phase tremor epochs.In that sense, asymmetric adaptivity might serve as a complex dynamical switching mechanism.In contrast, two anatomically intermingled populations (Fig. 1) activating the same muscle would display recurrent epochs of synchrony, interconnected by bursts of mass synchrony, potentially translating to pronounced tremor epochs, intersected by epochs of, e.g., less coherent twitches, lacking distinct tremor bursts.
Bursting, as it has been discussed in the literature, e.g., for individual neurons, is a phenomenon consisting of episodes of activity and inactivity. 34On the other hand, recurrent synchronization, as it has been introduced in this work, is a phenomenon arising from the interaction of oscillators in a complex dynamical network.It is characterized by bursting on the macroscopic level rather than on the level of individual microscopic dynamical units.
While this work is exemplified by using the Hodgkin-Huxley model to show the generality of our findings, the theory developed is not restricted to the field neuroscience.Moreover, we note that the model used in this work represents a phenomenological model for a complex dynamical system inspired by neuroscience rather than a realistic network model.As we show, recurrent synchronization is a result of the interplay of a coupling structure with the adaptation on a slower time scale.The "hidden" adaptation variables play the role of the slowly changing control parameters triggering the recurrent synchronization events.Recurrent synchronization induced by the Chaos ARTICLE scitation.org/journal/chainterplay of multiple timescales and complex dynamical networks is inherent in many dynamical networks modeling, e.g., social interactions or climate tipping dynamics. 70,71Bursting phenomena are also known dynamical states found outside the field of neuroscience. 72os θ dθ ω − c 1 sin θ − c 2 cos θ .(B1b) The integrals are analytically solvable, and we show it with the example of Eq. (B1a).We rewrite the integral term in the following form: Since the residues must lie within the 1-circle, we compute their absolute value; recall that A 2 = c 2 1 + c 2 2 , For ω > 0, the solution z + is located within the unit circle |z + | ≤ 1.We can, therefore, rewrite the first integral in Eq. (B2) in the following way: For the second integral in Eq. (B2), the same procedure yields We can now write Eq. (B2) as ω −

FIG. 2 .
FIG. 2. Recurrent synchronization in a network of Hodgkin-Huxley neurons.Panel a illustrates schematically the structure of the two-population network of Hodgkin-Huxley neurons with asymmetric spike-timing-dependent plasticity (shown in panel b).The colors of the network links refer to the adaptation rule implemented within and between the populations.Panel c shows the time series of the order parameter R (blue, left scale) displaying multiple cycles of bursting and phase-locked states.Panel d shows a zoom into the bursting episode.In panel e, the firing density is shown in blue for the whole network, in brown for population 1, and in green for population 2. Panels f and g display exemplary raster plots for the phase-locked and bursting episodes, respectively, with neurons of the first population (brown) and second population (green).For simplicity, we consider an all-to-all coupling structure on which the coupling weights evolve.The initial conditions for the coupling weights and the starting voltage were randomly chosen from the intervals [0, 0.5] and [−70 mV, 20 mV], respectively.The constant input currents are randomly chosen from the intervals [4.99 µA, 5.01 µA] and [12.99 µA, 13.01 µA] for populations 1 and 2. All other parameters are provided in Sec.II A.

FIG. 4 .FIG. 5 .
FIG. 4. Recurrent synchronization in a network of Hodgkin-Huxley neurons with heterogeneous constant currents.Time evolution of the order parameter for 200 Hodgkin-Huxley neurons after transient.The input currents for populations 1 and 2 are chosen from the intervals [4.5 µA, 5.5 µA] and [12.5 µA, 13.5 µA], respectively.The initial conditions for the coupling weights and the starting voltage were randomly chosen from the intervals [0, 0.5] and [−70 mV, 20 mV].

( 14 ) 6 ©FIG. 6 .
FIG.6.Recurrent synchronization in a system of two Hodgkin-Huxley neurons with asymmetric synaptic plasticity.Panel a shows a time series of the order parameter R, which exhibits recurrent synchronization.Panel b and c provide zooms into the time series of the order parameter for asynchronous spiking and at the transition from synchronous to asynchronous spiking, respectively.In panel d, the phase portrait of the reduced system in coupling variables (κ 1 , κ 2 ) is shown.Black lines show sample trajectories (flow lines) for the reduced system.The green curve is the projection of the whole system's trajectory onto the (κ 1 , κ 2 )-plane.The red line marks the boundary between the synchronous (light pink) and asynchronous (light blue) regimes.Panel e shows three different asymptotic states identified numerically in the parameter plane (τ 1 , γ ) with trajectories starting in (κ 1 , κ 2 ) = (0.15, 1).Asynchronous spiking is marked in blue, recurrent synchronization is shown in yellow and phase-locked spiking is depicted in red.The white areas indicate states that could not be categorized numerically.All parameters used for the simulation are given in Sec.II A.

FIG. 7 .
FIG. 7. Two Hodgkin-Huxley with identical update functions do not display recurrent synchronization.(a) and (b): reduced flow for two Hodgkin-Huxley neurons with identical update functions of the coupling weights; a displays the case of both coupling weights being updated by W 1 given by Eq. (12) of the main text; b shows the case of identical update functions W 2 given by Eq. (13) of the main text.The green curve is the projection of the trajectory of the whole system [starting in (κ 1 , κ 2 ) = (0.05, 1) in both cases] onto the (κ 1 , κ 2 )-plane and the red line marks the boundary between the synchronous (light pink) and asynchronous (gray blue) regions.(c)-(f): distributions of the timing difference between the spikes of neurons 1 and 2 (relative to the spiking times of neuron 1) for fixed coupling weights κ 2 = 0.1 and κ 1 = 0.05, 0.1, 0.2, and 0.25, which corresponds to the change from an asynchronous to a synchronous regime.

) Chaos 33 , 8 ©FIG. 8 .
FIG. 8. Recurrent synchronization in the phase oscillator model (4)-(7).(a), the bifurcation diagram in the (a, b)-parameter plane is shown.A and B denote the regions with stable recurrent synchronization, i.e., the stable limit cycle traversing the synchronous and asychronous domains in the (κ 1 , κ 2 ) plane.Region B corresponds to the coexistence of the bursting cycle and the stable equilibrium at κ 1 = κ 2 = 0. Panels (b) and (d) show the trajectories of the reduced system (19) in the (κ 1 , κ 2 )-plane that correspond to the time series in (c) and (e) and (f), respectively.The phase portrait d is complemented by an additional trajectory (blue) separating the basins of attraction of the recurrent synchronization limit cycle and the stable equilibrium at (κ 1 , κ 2 ) = (0, 0).The red line denotes the boundary between the asynchronous (gray blue) and the phase-locked (light pink) regions.Panels (c) (a = 0.5, b = 0.07) , (e), and (f) (a = 0.385, b = 0.125) show the time evolution of the order parameter R for recurrent synchronization (green) and the trajectory converging to the equilibrium at (0, 0) (dark red).Parameters: = 0.0001, ω = 0.1, α = 0.25π, β = −π/2.

FIG. 9 .
FIG. 9. Bifurcation diagram the reduced system(10).Bifurcation diagram in the parameter plane (a, b) for ω 0.1 and α = 0.25 π.The regions G1 and G2 (blue) mark areas where the equilibrium (κ 1 , κ 2 ) = (0, 0) is stable; C (red): a stable equilibrium on the critical manifold; A (orange): a stable recurrent synchronization limit cycle; E (green): a stable limit cycle on the critical manifold.The other regions mark areas of bistability of the equilibria described above.Also shown are analytically obtained bifurcation lines for the whole system (red) and the (0,0) equilibrium in the reduced system (blue), with SN denoting a Saddle-Node bifurcation, PF a Pitchfork bifurcation and H a Hopf bifurcation (dashed lines).The boundaries for regions A, B, E, and F are obtained numerically.The diagrams in the two columns on the right and the row in the bottom show exemplary trajectories for the reduced system in the (κ 1 , κ 2 ) plane for the regions denoted in the bifurcation diagram.Black lines mark stable trajectories, blue lines mark unstable trajectories, and green lines mark the limit cycle.The (0,0) equilibrium can either be stable (red and filled) or unstable (red and void).The blue areas denote the basins of attraction of the (0,0) equilibrium in the case of bistability.Green dots indicate stable fixed points on the critical manifold and the red line marks the boundary of the critical manifold.