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 fill 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 results 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

A multitude of real-world dynamical networks possess adaptive interactions.^{1–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–13} contributes to learning^{14} or anti-kindling^{15} effects. Another adaptation mechanism, structural plasticity, is responsible for a homeostatic regulation of electrical activity in the brain.^{16} Beyond neural networks, adaptive networks are used in communication systems^{17} and for modeling complex behavior in social, climate, or ecological systems.^{5,18,19}

In 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,21} So far, little is known about the dynamical effects induced by heterogeneous adaptation.^{22,23} The 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–27} where single neurons can have alternating periods of quiescence and fast spiking.

Other studies^{28–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.^{34} In 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.^{35} In 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.^{36}

In 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 re-emergence of episodes of synchronization and desynchronization. Effects such as converse symmetry breaking^{37,38} or asymmetry-induced 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.

## II. MODELS AND MEASURES

### A. Network of Hodgkin–Huxley neurons

^{40,41}

^{42}and of the following form:

^{43}and of the following form:

### B. Two phase oscillators

^{44,45}Taking into account the adaptation, the model of two adaptively coupled phase oscillators reads

^{46,47}

^{48}The parameters $a$ and $b$ scale the influence of the phase dynamics on the coupling weights. An additional asymmetry in the adaption rules for $ \kappa 1$ and $ \kappa 2$ is introduced by the parameter $\beta $.

### C. Macroscopic observables

^{44}Particularly, for Hodgkin–Huxley neurons considered in Sec. II A, each node possesses a sequence of spiking events ${ t i , 1,t i , 2,\u2026}$. The mapping used in this study is given by

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,52} These 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,55}

To 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 $\alpha $-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

^{56}

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

Using the explicit splitting of the time scales (smallness of $\u03f5$), we develop a further reduction in the mean-field equations (10)–(13). In particular, as the coupling weights $( \kappa 1, \kappa 2)$ evolve slowly, we consider them as parameters for the dynamics of the oscillators. Depending on the coupling weights, the coupled oscillator system $ x \u02d9 1= F 1( x 1, x 2, \kappa 1)$, $ x \u02d9 2= F 2( x 2, x 1, \kappa 2)$ shows either synchronized or desynchronized motion.

As the system (14) describes the (hidden) dynamics with respect to the slow time, it does not involve the small parameter $\u03f5$ 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 semi-analytically.

^{57}If the fast subsystem has one or more equilibria, then they satisfy the following equations:

## 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 $ \kappa 1= \kappa 12$ is updated accordingly to the STDP rule with the update function $ W 1$ and $ \kappa 2= \kappa 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$.

Importantly, the fast systems have to be computed once in order to find the distributions $\rho (\xi )$ for every grid point in the $ \kappa 1, \kappa 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 $( \kappa 1, \kappa 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 $ \tau 1$ (parameter of $ W 1$) and $\gamma $ (parameter of $ W 2$) is presented in Fig. 6(e). The parameter regions for which recurrent synchronization is observed are shown in yellow. The emergence of recurrent synchronization for different parameter regimes is clearly visible, with even two unconnected parametric regions. Therefore, the phenomenon of 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 $ \kappa 1$ and $ \kappa 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 $( \kappa 1, \kappa 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

The functions $ W ^ 1 , 2( \kappa 1, \kappa 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 $ \kappa 1$, $ \kappa 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 $\omega = \omega 1\u2212 \omega 2$ and $\alpha $. 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 $( \kappa 1, \kappa 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 $\beta $, see Fig. 12, and time-scale separation parameter $\u03f5$, see Fig. 11. In particular, the symmetric case $\beta =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 cycle 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 $( \kappa 1, \kappa 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 $\u03f5$, 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 color-coded 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 $( \kappa 1, \kappa 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,23}

Our 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. So far, the mechanism underlying the generation of Parkinsonian resting tremor remains an open question.^{63–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.^{66} Furthermore, causal data analysis suggests that this activity drives the muscular tremor, as opposed to being just driven by the sensory feedback from the periphery.^{67} In 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.^{65} According 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,68} However, 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,69} Based 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.^{34} On 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 interplay of multiple timescales and complex dynamical networks is inherent in many dynamical networks modeling, e.g., social interactions or climate tipping dynamics.^{70,71} Bursting phenomena are also known dynamical states found outside the field of neuroscience.^{72}

## ACKNOWLEDGMENTS

This work was supported by the German Research Foundation DFG, Project Nos. 411803875 and 440145547. P.A.T. gratefully acknowledges the funding support of this study by the Vaughn Bryson Research Fund and the John A. Blume Foundation.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Max Thiele:** Formal analysis (equal); Investigation (equal); Software (equal); Visualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Rico Berner:** Conceptualization (equal); Investigation (equal); Methodology (equal); Supervision (equal); Writing – review & editing (equal). **Peter A. Tass:** Investigation (equal); Writing – review & editing (equal). **Eckehard Schöll:** Investigation (equal); Project administration (equal); Supervision (equal); Writing – review & editing (equal). **Serhiy Yanchuk:** Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Project administration (equal); Supervision (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are openly available in GitHub, Ref. 73.

### APPENDIX A: ADIABATIC ELIMINATION FOR TWO PHASE OSCILLATORS

### APPENDIX B: AVERAGING OF THE TWO PHASE OSCILLATOR SYSTEM

### APPENDIX C: BIFURCATION ANALYSIS OF TWO COUPLED PHASE OSCILLATORS

## REFERENCES

*Neuronal Dynamics: From Single Neurons to Networks and Models of Cognition*

*Synchronization: A Universal Concept in Nonlinear Sciences*

*Chemical Oscillations, Waves and Turbulence*

*Averaging Methods in Nonlinear Dynamical Systems*