Chemical oscillators synchronized via an active oscillating medium: dynamics and phase approximation model

Chaos, Solitons and Fractals 145 (2021) 110809 Contents lists available at ScienceDirect Chaos, Solitons and Fractals Nonlinear Science, and Nonequilibrium and Complex Phenomena journal homepage: www.elsevier.com/locate/chaos Frontiers Chemical oscillators synchronized via an active oscillating medium: Dynamics and phase approximation model David García-Selfa a , b , Gourab Ghoshal c , Christian Bick d , Juan Pérez-Mercader e , f , Alberto P. Muñuzuri a , ∗ a CRETUS Institute, Group of Nonlinear Physics, Dept. of Physics, University of Santiago, Campus SurSantiago de Compostela, 15782, Spain b CESGA (Supercomputing Center of Galicia), Avda. de Vigo s/n, Santiago de Compostela, 15705, Spain c Department of Physics & Astronomy, University of Rochester, Rochester, New York, 14607 USA d Centre for Systems Dynamics and Control and Department of Mathematics, University of Exeter, Exeter, EX4 4QF, UK e Department of Earth and Planetary Sciences, Harvard University, Cambridge, 02138, Massachusetts, USA f The Santa Fe Institute, 1399 Hyde Park Road, Santa Fe, New Mexico, 87501, USA a r t i c l e i n f o Article history: Received 4 February 2021 Accepted 21 February 2021 Available online 6 March 2021 Keywords: Chemical oscillators Synchronization Bifurcation theory Phase approximation Mobbing state a b s t r a c t Different types of synchronization states are found when non-linear chemical oscillators are embedded into an active medium that interconnects the oscillators but also contributes to the system dynamics. Using different theoretical tools, we approach this problem in order to describe the transition between two such synchronized states. Bifurcation and continuation analysis provide a full description of the parameter space. Phase approximation modeling allows the calculation of the oscillator periods and the bifurcation point. ©2021 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license ( http://creativecommons.org/licenses/by-nc-nd/4.0/ ) 1. Introduction Chemical oscillatory behavior is evidence of complex, highly nonlinear dynamics, and is ubiquitous in nature [1,2] . In many cases, sub-units exhibiting oscillatory behavior couple together in large assemblies giving rise to a collective behavior of which a particularly important phenomenon is synchronization. Synchronization plays important roles in multiple biological and technical settings, for instance in the synchronized flashing of fireflies [3] , in cardiac pace-makers [4–6] , in yeast cells [7] , the firing of neurons [8] , in arrays of Josephson junctions [9] and semiconductor lasers [10] among numerous other examples. Given its ubiquity, the mechanism involved in synchronization has been the subject of rigorous study through both analytical considerations (phase models based on the Kuramoto-family of models [11,12] ) as well as experimental realizations including coupled electrochemical oscillators and reactors [13] . Of particular note are populations of catalyst-loaded oscillatory beads. A typical setup consists of a large number of beads in which the oscillatory Belousov–Zhabotinsky (BZ) reaction takes place [14] . These ∗Corresponding author. E-mail address: [email protected] (A.P. Muñuzuri). beads are immersed in a well-stirred active medium which acts as coupling between the population [15] , with oscillations being triggered by the contact of the beads with the medium. Several interesting dynamical behaviors have been reproduced in this setting including phase synchronization [16] , quorum sensing [14] and amplitude entrainment. Most studies carried out this far focus on the behavior of the beads while treating the active medium as just a mean to couple them, while ignoring its role as a potential oscillator in itself. Recently, however, a relatively new synchronization phenomenon was reported in numerical and experimental investigations, where the beads as well as the active medium are driven into a common high amplitude, low-frequency super-synchronized state of oscillations [17] . Interestingly, this occurred in the strong coupling limit, where previously the only reported state was that of oscillator death [14] . Beyond an experimental setting, this phenomenon has practical relevance, being similar to exotic states of synchronization as found in Interictal Epileptogenic Discharges, a known neuropathology [18,19] , and temperature mediated synchronization of the chirping of crickets [20] . Qualitative arguments and numerical analysis suggested the presence of higher harmonics [17] in the coupling function between the beads and the medium, although the precise forms were never presented. https://doi.org/10.1016/j.chaos.2021.110809 0960-0779/© 2021 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license ( http://creativecommons.org/licenses/by-nc-nd/4.0/ ) D. García-Selfa, G. Ghoshal, C. Bick et al. Chaos, Solitons and Fractals 145 (2021) 110809 In this paper, we fill this gap, by shedding light on the dynamical mechanisms behind this super-synchronized state of oscillations between the beads and the medium. In our theoretical approach, we consider a reduced system with two interacting oscillators, the collection of beads (that are synchronized a priori through standard coupling) and the active medium itself. The active medium is catalyst free and always coupled to all other oscillators. This reduced approach enables us to uncover the bifurcation structure for the system, across most of its known dynamical states. Reducing the dynamics to a set of phase equations, we calculate for the first time the period of oscillations for this exotic state [12,21] . In addition, we present the precise form for the higher-harmonics in the coupling function between the two oscillatory systems. The paper is organized as follows. In Section 2 we introduce the model equations with the needed simplifications and present the obtained state space diagram using bifurcation and continuation analysis (in this paper, we use the term space state diagram instead of phase diagram , as usually used in dynamical systems theory, in order to avoid confusions with the phase of the oscillators following [12] ). Finally, in Section 3 , we use the phase approximation model to study the transition between the synchronized and the super-synchronized states (mobbing states), calculating interaction functions, the periods of the oscillations in both of the states and the Fourier expansion of the interaction functions where we can see how the Fourier modes change in the transition between the synchronized and the super-synchronized states. The manuscript concludes with a section presenting the conclusions of this work. This manuscript points out the importance of an active medium as the means to couple oscillators. In fact, the active connecting medium introduces a great variety of non-trivial behaviors that cannot be described nor understood without its active dynamics. 2. Dynamic study: Bifurcation and continuation analysis We consider a system of n beads coupled chemical oscillators. Each oscillator is a resin bead loaded with the catalyst of the oscillatory BZ reaction. These beads are immersed in a surrounding solution containing all the chemicals of the BZ reaction except for the catalysts that it is in only present on the surface of the beads. The reactor is a continuously stirred tank [14,17] . This system can be described by the following set of differential equations. The dynamics of bead i ∈ { 1 , . . . , n beads } is described with the 3-variable Oregonator model [14,17] given by ⎧ ⎪ ⎨ ⎪ ⎩ ∂x i (t) ∂t = x i (t) ( 1 −x i (t) ) + y i ( t) ( q −x i (t) ) −K ex ( x i ( t) −x s (t) )  ∂y i (t) ∂t = 2 hz i (t) −y i (t) ( q + x i (t) ) −K ex ( y i (t) −y s (t) ) ∂z i (t) ∂t = x i (t) −z i (t) (1) where x i , y i and z i are the dimensionless variables representing the concentrations of activator, inhibitor, and catalyst, respectively for bead i . The quantities x s , y s are the dimensionless concentrations of activator and inhibitor in the surrounding solution (active medium), respectively. Note that the system is well-stirred and the concentration at any location of the surrounding medium is the same. The parameter K ex is the exchange rate constant between the beads and the surrounding solution. The parameters ,  , q and h are related to reaction rates and initial concentrations [17] . Since the surrounding solution (catalyst-free BZ reaction) interacts with all beads by exchanging activator and inhibitor in the reactor, it formally plays the role of coupling between the beads. For a well-stirred tank reactor, the dynamics of the concentrations of activator and inhibitor in the surrounding solution are given by ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ ∂x s (t) ∂t = x s (t) ( 1 −x s (t) ) + y s ( t) ( q −x s (t) ) +  V  n V s K ex n beads  i =1 ( x i (t) −x s (t) )  ∂y s (t) ∂t = −y s (t) ( q + x s (t) ) +  V  n V s K ex n beads  i =1 ( y i ( t) −y s (t) ) (2) where the parameter  V  n represents the average volume of the beads and the parameter V s is the total volume of the surrounding solution. Note that although the surrounding solution does not contain catalyst by itself, it does contain the beads that have the catalyst incorporated. Thus, the surrounding solution under these circunstances can potentially exhibit oscillations. In the following, we will focus in understanding the transition between the synchronized state to the mobbing state, i.e., the transition from the state characterized for all the chemical oscillators oscillating synchronous to the state in which all the beads oscillate in synchrony and with the surrounding medium also oscillating with the same amplitude and frequency. Given that the transitions of interest require all beads to be in a synchronized state, we consider that all oscillators, excluding the active medium, are identical and, thus, we can consider the synchronization manifold where the state of all oscillators is equal. Specifically, we assume that x i = x b , y i = y b , z i = z b for all i ∈ { 1 , . . . , n beads } (the index bdenotes a bead). Thus, on the synchronization manifold, the model (1),(2) reduce to the five-dimensional system ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂x b (t) ∂t = x b (t) ( 1 −x b (t) ) + y b ( t) ( q −x b (t) ) −K ex ( x b (t) −x s (t) )  ∂y b (t) ∂t = 2 hz b (t) −y b (t) ( q + x b (t) ) −K ex ( y b (t) −y s (t) ) ∂z b (t) ∂t = x b (t) −z b (t) ∂x s (t) ∂t = x s (t) ( 1 −x s (t) ) + y s ( t) ( q −x s (t) ) + ρK ex ( x b (t) −x s (t) )  ∂y s (t) ∂t = −y s (t) ( q + x s (t) ) + ρK ex ( y b ( t) −y s ( t) ) (3) with ρ= n beads  V  n V s is the density of the system. State space diagram For the five-dimensional simplified system, we obtain the state space diagram using continuation and bifurcation analysis of dynamical systems software ( Matcont [22] and AUTO [23] ). In Fig. 1 we show the state space diagrams obtained with the model in Eq. (3) . Fig. 1 a displays all the observed behaviors as a function of the exchange rate constant between the beads and the surrounding medium ( K ex ) and the density of beads ( ρ). As a first observation it is noteworthy that the same behaviors observed both experimentally and numerically are also observed with the same distribution on the state space diagram [17] . The same bifurcation diagram is plotted in Fig. 1 b in a 3D perspective where the vertical axis corresponds with the value of the variable x b for the beads. This new representation reveals the details of the different bifurcations involved in the transitions analyzed. In both representations we observe, the generalized Hopf bifurcation point ( GH) separating the two branches of supercritical Hopf bifurcation ( H −), where the first Lyapunov coefficient is negative, and subcritical Hopf bifurcation ( H + ), where the first Lyapunov coefficient is positive, and the saddle-node bifurcation of periodic orbits ( LP Ccurve), where the system has a unique non-hyperbolic limit cycle with the nontrivial Floquet multiplier +1 . On the other hand, coincident with the GH point, we have a cusp point of cycles ( C P C ) as this bifurcation sepa- 2 D. García-Selfa, G. Ghoshal, C. Bick et al. Chaos, Solitons and Fractals 145 (2021) 110809 Fig. 1. State space diagram. (a) Different behaviors observed in the system when K ex and ρare varied. (b) Same state space diagram but in a three-dimensional perspective with the value os x b at the stationary is plotted in the vertical axis. Fig. 2. (a) Limit cycles corresponding to the transitions from equilibrium point A (oscillation death) in point A to super-synchronized state in point B passing through the synchronized state (between the supercritical Hopf bifurcation, H, and the saddle-node bifurcations of periodic orbits, LPC). Dashed green line A-B corresponding to that shown in Fig. 1 a. (b) Frequencies corresponding to these transitions. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) rates the supercritical behavior from the subcritical one. This simplified model captures the dynamics of the system that was described in Figure 2-d of [17] . Note that, since all the beads are identical by construction in our system, the non-synchronization phase shown in [17] does not appear. In Fig. 2 a we observe periodic orbits in the transitions from equilibrium (this state is the equivalent to the oscillations death in the complete ( 3 ×n beads + 2 )-dimensional model) in point A to the super-synchronized state in point B, passing through the supercritical Hopf Hand the saddle-node bifurcations of periodic orbits LP C. Note that the rapid change of the limit cycle is a Canard explosion that arises even in the simple Oregonator model [24–26] . We also observe a drastic decrease in frequency in the Canard explosion in Fig. 2 b. Thus, the Canard explosion observed in a simple Oregonator can also be observed in the synchronized population of Oregonators coupled via the active medium. 3. Phase approximation model In the following, we use phase approximation in order to obtain the periods of the beads in the synchronized and supersynchronized phases. Given an oscillating system, its state is described by its position along its limit cycle (its phase). If we have two uncoupled oscillators, their phases lie on a torus. If the two oscillators have stable limit cycles and they are weakly coupled, the torus persists and we can describe their state by their phases [12,21] . Theoretical model Consider two oscillators coupled through an interaction function G i such that the dynamics are given by d X 1 (t) dt = F ( X 1 ) + K G 1 ( X 1 , X 2 ) d X 2 (t) dt = F ( X 2 ) + K G 2 ( X 1 , X 2 ) (4) If the coupling strength Kis small, the dynamics can be reduced to a phase description, that is, the state of each oscillator is determined by a phase variable θi , i ∈ { 1 , 2 } on the circle which evolves according to dθ1 (t) dt = ω 1 + KH 1 (θ2 −θ1 ) dθ2 (t) dt = ω 2 + KH 2 (θ1 −θ2 ) (5) where ω i are the intrinsic frequencies of the oscillators and H i (φ) are the phase interaction functions (that depend only on the phase difference φ= θ2 −θ1 ). The phase interaction function is computed by averaging %beginwidetext H i (φ) = 1 T  T 0 Z (t) ·G i ( X 0 (t + φ) , X 0 (t)) dt = 1 2 π 2 π 0 Z (ϕ) ·G i ( X 0 (ϕ + φ) , X 0 (ϕ )) dϕ (6) where X 0 (t) is stable limit cycle, Z (t) is the adjoin or phase response curve (PRC)—the phase shift function obtained when the system that lies on its limit cycle is infinitesimally perturbed—and ·denotes the scalar product of vectors. Both Z (θ)and H i (φ)can be 3 D. García-Selfa, G. Ghoshal, C. Bick et al. Chaos, Solitons and Fractals 145 (2021) 110809 obtained numerically using XPPAUTO [23] . For a theoretical derivation of the PRC by means of the adjoint method and the average method for calculation of the interaction function, see [12] . Application to our problem The problem considered, consists on a bead (all the initial beads, with density ρ, are considered to be completely synchronized-) uncoupled to the rest of the system (the surrounding solution). Note that we are assuming that all the beads are already in a synchronized state and, thus, can be represented by one single set of equations as in previous section. This system will be numerically solved using the above mentioned software. The equations describing two identical copies of the system are ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂x b i (t) ∂t = 1 x b i (t) 1 −x b i (t) + y b i ( t) q −x b i (t)  ∂y b i (t) ∂t = 1  2 hz b i (t) −y b i (t) q + x b i (t)  ∂z b i (t) ∂t = x b i (t) −z b i (t) ∂x s i (t) ∂t = 1 ( x s i (t) ( 1 −x s i (t) ) + y s i ( t) ( q −x s i (t) ) + ρK ex x b i (t) −x s i (t)  ∂y s i (t) ∂t = 1  −y s i (t) ( q + x s i (t) ) + ρK ex y b i ( t) −y s i ( t)  (7) with i = 1 , 2 . For the set of parameters in the oscillatory regime as described above, the beads and the surrounding solution of each copy oscillate with a natural frequency ω 0 . We now couple the two systems by coupling the beads of one system to the solution of the second system. Specifically, we couple the bead of system i = 1 with the surrounding solution of system i = 2 using G 1 = ⎡ ⎢ ⎢ ⎢ ⎣ −1 x b 1 −x s 2  −1  y b 1 −y s 2  0 0 0 ⎤ ⎥ ⎥ ⎥ ⎦ (8) and with a coupling strength K = K ex 2 ( 1+ ρ) . This allows us to derive a phase description if the coupling strength Kis small. For simplicity, we absorb the coupling strength into the coupling function H i so that the phase of oscillator 1 evolves according to dθ1 (t) dt = ω 0 + H 1 (θ2 −θ1 ) . (9) Now, we can calculate the periods of the oscillators in both synchronized and super-synchronized states and compare the full nonlinear model and the phase approximation. For the phase dynamics we know that the beads and the surrounding solution oscillate in phase, φ0 = 0 = 2 π, with frequency and period T = 2 π  so that dθ1 (t) dt = = 2 π T = ω 0 + H 1 (0) (10) For the full nonlinear system, the model chemical parameters were set q = 0 . 002 , = 0 . 01 ,  = 0 . 015 , h = 0 . 70 (11) The control parameters considered for the analysis are the chemical exchange rate K ex and the density of beads in the medium ρ. The components of the PRC Z (θ) ) for one bead with a density ρ= 1 . 2 are shown in Fig. 3 . Interaction functions ( H 1 (φ) ) in the synchronized and supersynchronized states and density ρ= 1 . 2 are shown in Fig. 4 . The phase approximation yields a good description of the period of the collective oscillation as parameters are varied. The calculated periods for each value of K ex obtained via numerical integration and by using the phase approximation model (10) with Fig. 3. Phase response curves (PRC) for one bead with ρ= 1 . 2 . Chemical parameters set used: q = 0 . 002 , = 0 . 01 ,  = 0 . 015 , h = 0 . 70 . Fig. 4. Two cases of H 1 (φ) with density ρ= 1 . 2 corresponding to (a) synchronized state with K ex = 0 . 033 and (b) super-synchronized state with K ex = 0 . 077 . Note the remarkable difference in scale on the y -axis of both figures. density ρ= 1 . 2 are shown in Fig. 5 . A transition from synchronized to super-synchronized phase is seen, with a discontinuity in the period and the interaction function caused by the bifurcation described between limit cycles of different nature. 4 D. García-Selfa, G. Ghoshal, C. Bick et al. Chaos, Solitons and Fractals 145 (2021) 110809 Fig. 5. Periods of the oscillators in function of K ex with density ρ= 1 . 2 . The discontinuity in period shows the transition from synchronized to super-synchronized states. Fig. 6. Fourier coefficients (sine-cosine serie). Due to the different nature of the limit cycles in the states of synchronization and super-synchronization, we have different modes in the function of interaction in these states. Fourier expansion The qualitative change in the oscillation in the transition between the synchronized and super-synchronized states can also be seen in the change of interaction function of the phase reduction. Indeed, the changes in the phase interaction function are an indicator of the underlying bifurcations [27] . To illustrate this effect in the chemical oscillator system, we expand the interaction function in Fourier series, and to understand how the Fourier modes change as the system parameters are varied. We can expand the interaction function into a sine-cosine Fourier series (in the supplementary information a exponential Fourier series expansion is presented showing equivalent results), H 1 (φ) = a 0 2 + ∞  k =1 [ a k cos ( kφ) + b k sin ( kφ) ] (12) In Fig. 6 , we show the coefficients obtained for each value of K ex ( H 1 (φ) calculated with the method referred above for each value of K ex ). We can see that in the synchronized state: H 1 (φ) ∼a 0 2 + a 1 cos ( φ) + b 1 sin ( φ) , (13) and the values of the coefficients hardly vary, except in the vicinity of the transition between synchronized and super-synchronized states. On the other hand, in the super-synchronized state we have higher harmonics and the value of the coefficients varies considerably for a wide range of values of K ex : H 1 (φ) ∼a 0 2 + a 1 cos ( φ) + b 1 sin ( φ) + a 2 cos ( 2 φ) + b 2 sin ( 2 φ) + a 3 cos ( 3 φ) + b 3 sin ( 3 φ) (14) 4. Conclusions In this manuscript, we considered the problem of synchronization between oscillators embedded into an active medium. This type of system has shown to exhibit more than one state of synchronization. The chosen system is constituted by a set of chemical oscillators immersed into a chemical solution that provides the physical medium to interact but also adds dynamics to the total system. This system has been shown to synergetically produce a different synchronization state (supersynchronization) that was not accessible for each of the two main components of the problem (external medium or the beads). Using continuation and bifurcation analysis, we reconstructed the experimental and numerical parameter space previously reported in [17] . Three different states of synchronization are found; oscillation death, synchronization and mobbing state (supersynchronization) and the transitions between each other analyzed from a bifurcation analysis point of view. Note that in the fivedimensional simplified model used, the oscillation death correspondes with a steady state of the system, synchronization is a normal oscillatory behavior and supersynchronized state is demonstrated as a different state of oscillation. In order to calculate the periods exhibited by the oscillators in each state, a phase approximation model was considered reproducing with good accuracy the previous results reported both in experiments and in numerical simulations as well as the discontinuity that signals the transition between synchronization to mobbing state. Finally, we have proved that the discontinuity in the periods and the interaction function in the transition from synchronized to super-synchronized states is caused by the bifurcation described between limit cycles of different nature. This was possible to understand considering that the active medium is an oscillator although with a different nature (the medium does not have a catalyst per se although the catalyst is included into the beads that are immersed in the medium and, thus, needs the activity of the beads to oscillate). The results of this analysis can extrapolate to different systems as far as the connective medium plays an active role in the dynamics of the system showing the generality of the phenomenon described. This type of system is found in different fields in Nature including neuronal processes involving the glia [28] , as glial cells and neurons have ionic channels that allow them to oscillate but only neurons possess synaptic connections and have the ability to oscillate by themselves. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. 5 D. García-Selfa, G. Ghoshal, C. Bick et al. Chaos, Solitons and Fractals 145 (2021) 110809 Acknowledgments We gratefully acknowledge financial support by the Spanish Ministerio de Economía y Competitividad and European Regional Development Fund under contract RTI2018-097063-B-I00 AEI/FEDER , UE, and by Xunta de Galicia under Research Grant No. 2018-PG082 . APM and DGS are part of the CRETUS Strategic Partnership (AGRUP2015/02). All these programs are co-funded by FEDER (UE). JP-M thanks Repsol S. A. for supporting this research. Funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. References [1] Strogatz S . Sync: the emerging science of spontaneous order. Hyperion Press; 2003 . ISBN 0786 86 844 9 [2] Murray JD. Mathematical biology. 2nd. Springer-Verlag Berlin Heidelberg; 1993. doi: 10.1007/978- 3- 662- 08542- 4 . [3] Buck J. Synchronous rhythmic flashing of fireflies. II.. Q Rev Biol 1988;63(3):265–89 . http://www.jstor.org/stable/2830425 . [4] Peskin CS . Mathematical aspects of heart physiology. New York: Courant Institute of Mathematical Sciences, New York University; 1975 . [5] Torre V. A theory of synchronization of heart pace-maker cells. J Theor Biol 1976;61(1):55–71. doi: 10.1016/0022-5193(76)90104-1 . [6] Mirollo RE, Strogatz SH. Synchronization of pulse-coupled biological oscillators. SIAM J Appl Math 1990;50(6):1645–62. doi: 10.1137/0150098 . [7] Ghosh A, Chance B, Pye E. Metabolic coupling and synchronization of NADH oscillations in yeast cell populations. Arch Biochem Biophys 1971;145(1):319– 31. doi: 10.1016/0 0 03-9861(71)90 042-7 . [8] Izhikevich EM . Dynamical systems in neuroscience: the geometry of excitability and bursting. Cambridge, Mass: MIT Press; 2007 . [9] Wiesenfeld K, Colet P, Strogatz SH. Frequency locking in Josephson arrays: connection with the kuramoto model. Physical Review E - Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics 1998;57(2):1563–9. doi: 10. 1103/PhysRevE.57.1563 . [10] Hohl A, Gavrielides A, Erneux T. Localized synchronization in two coupled nonidentical semiconductor lasers. Phys Rev Lett 1997;78(25):4745–8. doi: 10. 1103/PhysRevLett.78.4745 . [11] Acebrón JA, Bonilla LL, Pérez Vicente CJ, Ritort F, Spigler R. The kuramoto model: a simple paradigm for synchronization phenomena. Rev Mod Phys 2005;77(1):137–85. doi: 10.1103/RevModPhys.77.137 . [12] Nakao H. Phase reduction approach to synchronisation of nonlinear oscillators. Contemp Phys 2016;57(2):188–214. doi: 10.1080/00107514.2015.1094987 . [13] Kiss IZ, Zhai Y, Hudson JL. Predicting mutual entrainment of oscillators with experiment-Based phase models. Phys Rev Lett 2005;94(24):248301. doi: 10. 1103/PhysRevLett.94.248301 . [14] Taylor AF, Tinsley MR, Wang F, Huang Z, Showalter K. Dynamical quorum sensing and synchronization in large populations of chemical oscillators. Science 2009;323(5914) 614LP–617. doi: 10.1126/science.1166253 . [15] Tinsley M, Taylor A, Huang Z, Wang F, Showalter K. Dynamical quorum sensing and synchronization in collections of excitable and oscillatory catalytic particles. Physica D 2010;239(11):785–90. doi: 10.1016/J.PHYSD.20 09.08.0 01 . [16] Taylor AF, Tinsley MR, Wang F, Showalter K. Phase clusters in large populations of chemical oscillators. Angewandte Chemie 2011;123(43):10343–6. doi: 10.1002/ange.201008248 . [17] Ghoshal G, Muñuzuri AP, Pérez-Mercader J. Emergence of a super-synchronized mobbing state in a large population of coupled chemical oscillators. Sci Rep 2016;6:19186. doi: 10.1038/srep19186 . [18] Prince DA , Connors BW . Mechanisms of interictal epileptogenesis.. Adv Neurol 1986;44:275–99 . [19] de Curtis M, Avanzini G. Interictal spikes in focal epileptogenesis. Prog Neurobiol 2001;63(5):541–67. doi: 10.1016/S0301-0 082(0 0)0 0 026-5 . [20] Walker TJ. Acoustic synchrony: two mechanisms in the snowy tree cricket. Science 1969;166(3907):891–4. doi: 10.1126/science.166.3907.891 . [21] Pietras B, Daffertshof er A. Network dynamics of coupled oscillators and phase reduction techniques. Phys Rep 2019;819:1–105. doi: 10.1016/j.physrep.2019.06. 001 . [22] Dhooge A, Govaerts W, Kuznetsov YA, Meijer HGE, Sautois B. New features of the software MatCont for bifurcation analysis of dynamical systems. Math Comput Model Dyn Syst 2008;14(2):147–75. doi: 10.1080/13873950701742754 . [23] Ermentrout B . Simulating, analyzing, and ani-mating dynamical systems: a guide toXPPAUT for researchers and students. Philadelphia: SIAM; 2002 . ISBN 0-89871-506-7. [24] Bo Peng, Vilmos G, Showalter K. False bifurcations in chemical systems: canards. Philosophical Transactions of the Royal Society of London Series A: Physical and Engineering Sciences 1991;337(1646):275–89. doi: 10.1098/rsta. 1991.0123 . [25] Brøns M, Bar-Eli K. Canard explosion and excitation in a model of the Belousov-Zhabotinsky reaction. J Phys Chem 1991;95(22):8706–13. doi: 10. 1021/j100175a053 . [26] Krupa M, Szmolyan P. Relaxation oscillation and canard explosion. J Differ Equ 2001;174(2):312–68. doi: 10.1006/jdeq.2000.3929 . [27] Hesse J, Schleimer J-H, Schreiber S. Qualitative changes in phase-response curve and synchronization at the saddle-node-loop bifurcation. Physical Review E 2017;95(5):52203. doi: 10.1103/PhysRevE.95.052203 . [28] Alvarez-Maubecin V, García-Hernández F, Williams JT, Van Bockstaele EJ. Functional coupling between neurons and glia. The Journal of Neuroscience 2018;20(11):4091–8. doi: 10.1523/jneurosci.20- 11- 04091.20 0 0 . 6