Chemical oscillators synchronized via an active oscillating medium: dynamics and phase approximation model
David García-Selfa, Gourab Ghoshal, Christian Bick, Juan Pérez-Mercader, Alberto P. Muñuzuri
aa r X i v : . [ n li n . AO ] F e b Chemical oscillators synchronized via an active oscillating medium: dynamics andphase approximation model
David Garc´ıa-Selfa,
Gourab Ghoshal, Christian Bick, Juan P´erez-Mercader, and Alberto P. Mu˜nuzuri CRETUS Institute. Group of Nonlinear Physics. Dept. of Physics.Universidade de Santiago de Compostela. 15782 Santiago de Compostela, Spain CESGA (Supercomputing Center of Galicia). 15705 Santiago de Compostela, Spain Department of Physics & Astronomy, University of Rochester, Rochester, New York 14607, USA Centre for Systems Dynamics and Control and Department of Mathematics, University of Exeter, Exeter EX4 4QF, UK Department of Earth and Planetary Sciences. Harvard University, Cambridge, MA 02138, USA The Santa Fe Institute, 1399 Hyde Park Road, Santa Fe, NM 87501, USA
Different types of synchronization states are found when non-linear chemical oscillators are em-bedded into an active medium that interconnects the oscillators but also contributes to the systemdynamics. Using different theoretical tools, we approach this problem in order to describe the tran-sition between two such synchronized states. Bifurcation and continuation analysis provide a fulldescription of the parameter space. Phase approximation modeling allows the calculation of theoscillator periods and the bifurcation point.
I. 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 collectivebehavior of which a particularly important phenomenon is synchronization. Synchronization plays important rolesin 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 semiconductorlasers [10] among numerous other examples.Given its ubiquity, the mechanism involved in synchronization has been the subject of rigorous study through bothanalytical considerations (phase models based on the Kuramoto-family of models [11, 12]) as well as experimentalrealizations including coupled electrochemical oscillators and reactors [13]. Of particular note are populations ofcatalyst-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 beads are immersed in a well-stirred active medium which acts ascoupling 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 amean to couple them, while ignoring its role as a potential oscillator in itself. Recently, however, a relatively newsynchronization phenomenon was reported in numerical and experimental investigations, where the beads as well as theactive 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 oscillatordeath [14]. Beyond an experimental setting, this phenomenon has practical relevance, being similar to exotic states ofsynchronization as found in Interictal Epileptogenic Discharges, a known neuro-pathology [18, 19], and temperaturemediated synchronization of the chirping of crickets [20]. Qualitative arguments and numerical analysis suggested thepresence of higher harmonics [17] in the coupling function between the beads and the medium, although the preciseforms were never presented.In this paper, we fill this gap, by shedding light on the dynamical mechanisms behind this super-synchronizedstate of oscillations between the beads and the medium. In our theoretical approach, we consider a reduced systemwith 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. Thisreduced approach enables us to uncover the bifurcation structure for the system, across most of its known dynamicalstates. Reducing the dynamics to a set of phase equations, we calculate for the first time the period of oscillations forthis exotic state [12, 21]. In addition, we present the precise form for the higher-harmonics in the coupling functionbetween the two oscillatory systems.The paper is organized as follows. In Section II we introduce the model equations with the needed simplifications andpresent 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 confusionswith the phase of the oscillators following [12]). Finally, in Section III, we use the phase approximation model to studythe transition between the synchronized and the super-synchronized states (mobbing states), calculating interactionfunctions, the periods of the oscillations in both of the states and the Fourier expansion of the interaction functionswhere we can see how the Fourier modes change in the transition between the synchronized and the super-synchronizedstates. The manuscript concludes with a section presenting the conclusions of this work.This manuscript really 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 neitherunderstood without its active dynamics.
II. 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 catalystof the oscillatory BZ reaction. These beads are immersed in a surrounding solution containing all the chemicals of theBZ reaction except for the catalysts that it is in only present on the surface of the beads. The reactor is a continuouslystirred tank [14, 17].This system can be described by the following set of differential equations. The dynamics of bead i ∈ { , ..., 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 thesurrounding solution (active medium), respectively. Note that the system is well-stirred and the concentration at anylocation of the surrounding medium is supposed to be the same. The parameter K ex is the exchange rate constantbetween the beads and the surrounding solution. The parameters ǫ , ǫ ′ , q and h are related to reaction rates and initialconcentrations [17].Since the surrounding solution (catalyst-free BZ reaction) interacts with all beads by exchanging activator andinhibitor in the reactor, it formally plays the role of coupling between the beads. For a well-stirred tank reactor, thedynamics 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 )) + h V i n V s K ex n beads X i =1 ( x i ( t ) − x s ( t )) ǫ ′ ∂y s ( t ) ∂t = − y s ( t ) ( q + x s ( t )) + h V i n V s K ex n beads X i =1 ( y i ( t ) − y s ( t )) (2)where the parameter h V i n represents the average volume of the beads and the parameter V s is the total volume ofthe surrounding solution. Note that although the surrounding solution does not contain catalyst by itself, it doescontain the beads that have the catalyst incorporated. Thus, the surrounding solution under these circunstances canpotentially 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 inwhich all the beads oscillate in synchrony and with the surrounding medium also oscillating with the same amplitudeand frequency.As the transitions we are interested in involve that all the beads are already in a synchronized state, we considerthat all oscillators, excluding the active medium, are identical and, thus, we can consider the synchronization manifoldwhere 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 ∈ { , . . . , n beads } (the index b denotes a bead). Thus, on the synchronization manifold, the model (1),(2) reduce to the five-dimensionalsystem ǫ ∂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 h V i 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 bifurcationanalysis of dynamical systems software (
Matcont [22] and
AUTO [23]). In
Figure 1 we show the state space diagramsobtained with the model in Eq. (3).
Figure 1a displays all the observed behaviors as a function of the exchange rateconstant between the beads and the surrounding medium ( K ex ) and the density of beads ( ρ ). As a first observation itis noteworthy the fact that the same behaviors observed both experimentally and numerically are also observed withthe same distribution on the state space diagram [17]. The same bifurcation diagram is plotted in Figure 1b in a 3Dperspective where the vertical axis corresponds with the value of the variable x b for the beads. This new representationunveils 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 − ), wherethe first Lyapunov coefficient is negative, and subcritical Hopf bifurcation ( H + ), where the first Lyapunov coefficientis positive, and the saddle-node bifurcation of periodic orbits ( LP C curve), 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 ( CP C ) as this bifurcation separates the supercritical behavior from the subcriticalone. 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 notappear. In
Figure 2a we can see the periodic orbits in the transitions from equilibrium (this state is the equivalent
GH CPCLPCH - H + LPCH- H+ AB (a) (b) FIG. 1. State space diagram. (a) Different behaviors observed in the system when K ex and ρ are varied. (b) Same state spacediagram but in a three-dimensional perspective with the value os x b at the stationary is plotted in the vertical axis. to the oscillations death in the complete (3 × n beads + 2)-dimensional model) in point A to super-synchronized statein point B passing through the synchronized state (between the supercritical Hopf bifurcation H and the saddle-node bifurcations of periodic orbits LP C ). Note that the rapid change of the limit cycle is a Canard explosion that x b H z b LPCLPCLPCLPC LPCLPC
A B (a) (b)
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 thesaddle-node bifurcations of periodic orbits, LP C ). Dashed green line A-B corresponding to that shown in
Figure 1a . (b) Frequencies corresponding to these transitions. arises even in the simple Oregonator model [24–26]. We can also see the drastic decrease in frequency in the Canardexplosion in
Figure 2b . Thus, the Canard explosion observed in a simple Oregonator can also be observed in thesynchronized population of Oregonators coupled via the active medium.
III. PHASE APPROXIMATION MODEL
In the following, we use phase approximation in order to obtain the periods of the beads in the synchronized andsuper-synchronized phases. Given an oscillating system, its state is described by its position along its limit cycle (itsphase). If we have two uncoupled oscillators, their phases lie on a torus. If the two oscillators have stable limit cyclesand 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 ( t ) dt = F ( X ) + K G ( X , X ) d X ( t ) dt = F ( X ) + K G ( X , X ) (4)If the coupling strength K is small, the dynamics can be reduced to a phase description, that is, the state of eachoscillator is determined by a phase variable θ i , i ∈ { , } on the circle which evolves according to dθ ( t ) dt = ω + KH ( θ − θ ) dθ ( t ) dt = ω + KH ( θ − θ ) (5)where ω i are the intrinsic frequencies of the oscillators and H i ( φ ) are the phase interaction functions (that dependonly on the phase difference φ = θ − θ ). The phase interaction function is computed by averaging H i ( φ ) = 1 T Z T Z ( t ) · G i ( X ( t + φ ) , X ( t )) dt = 12 π Z π Z ( ϕ ) · G i ( X ( ϕ + φ ) , X ( ϕ )) dϕ (6)where X ( t ) is stable limit cycle, Z ( t ) is the adjoin or phase response curve (PRC)—the phase shift function obtainedwhen 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 obtained numerically using XPPAUTO [23]. For a theoretical derivation of the PRC bymeans 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 completelysynchronized among them) uncoupled to the rest of the system (the surrounding solution). Note that we are assumingthat all the beads are already in a synchronized state and, thus, can be represented by one single set of equationsas in previous section. This system will be numerically solved using the above mentioned software. The equationsdescribing 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 surroundingsolution of each copy oscillate with a natural frequency ω .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 = − ǫ ( x b − x s ) − ǫ ′ ( y b − y s )000 (8)and with a coupling strength K = K ex ρ ) . This now allows to derive a phase description if the coupling strength K issmall. For simplicity, we absorb the coupling strength into the coupling function H i so that the phase of oscillator 1evolves according to dθ ( t ) dt = ω + H ( θ − θ ) . (9)Now, we can calculate the periods of the oscillators in both synchronized and super-synchronized states and comparethe full nonlinear model and the phase approximation. For the phase dynamics we know that the beads and thesurrounding solution oscillate in phase, φ = 0 = 2 π , with frequency Ω and period T = π Ω so that dθ ( t ) dt = Ω = 2 πT = ω + H (0) (10)For the full nonlinear system, the model chemical parameters were set q = 0 . , ǫ = 0 . , ǫ ′ = 0 . , h = 0 .
70 (11)The control parameters considered for the analysis are the chemical exchange rate K ex and the density of beads inthe medium ρ . The components of the PRC Z ( θ )) for one bead with a density ρ = 1 . Figure 3 .Interaction functions ( H ( φ )) in the synchronized and super-synchronized states and density ρ = 1 . Figure 4 . /2 3 /2 20369 Z x b () Z y b () (a)(b) FIG. 3. Phase response curves (PRC) for one bead with ρ = 1 .
2. Chemical parameters set used: q = 0 .
002 , ǫ = 0 .
01 , ǫ ′ = 0 . h = 0 . 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 approximationmodel (10) with density ρ = 1 . Figure 5 . A transition from synchronized to super-synchronized phaseis seen, with a discontinuity in the period and the interaction function caused by the bifurcation described betweenlimit cycles of different nature.
Fourier expansion
The qualitative change in the oscillation in the transition between the synchronized and super-synchronized statescan also be seen in the change of interaction function of the phase reduction. Indeed, the changes in the phaseinteraction function are an indicator of the underlying bifurcations [27]. To illustrate this effect in the chemicaloscillator system, we expand the interaction function in Fourier series, and to understand how the Fourier modes /2 3 /2 2-4-2024 H () K ex =0.077 -3.9-3.8-3.7-3.6-3.5 H () K ex =0.033 (a)(b) FIG. 4. Two cases of H ( φ ) with density ρ = 1 . K ex = 0 .
033 and (b) super-synchronized state with K ex = 0 . y -axis of both figures. change as the system parameters are varied.We can expand the interaction function into a sine-cosine Fourier series (in the supplementary information aexponential Fourier series expansion is presented showing equivalent results), H ( φ ) = a ∞ X k =1 [ a k cos ( kφ ) + b k sin ( kφ )] (12)In Figure 6 , we show the coefficients obtained for each value of K ex ( H ( φ ) calculated with the method referredabove for each value of K ex ). We can see that in the synchronized state: H ( φ ) ∼ a a cos ( φ ) + b sin ( φ ) , (13)and the values of the coefficients hardly vary, except in the vicinity of the transition between synchronized andsuper-synchronized states.On the other hand, in the super-synchronized state we have higher harmonics and the value of the coefficients variesconsiderably for a wide range of values of K ex : H ( φ ) ∼ a a cos ( φ ) + b sin ( φ ) + a cos (2 φ ) + b sin (2 φ ) + a cos (3 φ ) + b sin (3 φ ) (14) K ex T ( t. u . ) =1.2 Phase approximationNumerical integration
FIG. 5. Periods of the oscillators in function of K ex with density ρ = 1 .
2. The discontinuity in period shows the transitionfrom synchronized to super-synchronized states.
IV. CONCLUSIONS
In this manuscript, we considered the problem of synchronization between oscillators embedded into an activemedium. This type of system has shown to exhibit more than one state of synchronization. The chosen system isconstituted by a set of chemical oscillators immersed into a chemical solution that provides the physical medium tointeract but also adds dynamics to the total system. This system has been shown to synergetically produce a differentsynchronization state (supersynchronization) that was not accessible for each of the two main components of theproblem (external medium or the beads).Using continuation and bifurcation analysis, we reconstructed the experimental and numerical parameter spacepreviously reported in [17]. Three different states of synchronization are found; oscillation death, synchronization andmobbing state (super-synchronization) and the transitions between each other analyzed from a bifurcation analysispoint of view. Note that in the five-dimensional simplified model used, the oscillation death correspondes with a steadystate of the system, synchronization is a normal oscillatory behavior and supersynchronized state is demonstrated asa different state of oscillation.In order to calculate the periods exhibited by the oscillators in each state, a phase approximation model wasconsidered reproducing with good accuracy the previous results reported both in experiments and in numericalsimulations 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 fromsynchronized to super-synchronized states is caused by the bifurcation described between limit cycles of differentnature. This was possible to understand considering that the active medium is an oscillator although with a differentnature (the medium does not have a catalyst per se although the catalyst is included into the beads that are immersedin 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 activerole in the dynamics of the system showing the generality of the phenomenon described. This type of system is foundin different fields in Nature including neuronal processes involving the glia [28], as glial cells and neurons have ionicchannels that allow them to oscillate but only neurons possess synaptic connections and have the ability to oscillateby themselves. K e(cid:0) -80 -(cid:1)(cid:2) K (cid:23)(cid:24) -30 (cid:25)(cid:26)(cid:27) -10010 b (cid:28) (cid:29)(cid:30)(cid:31) !" FIG. 6. Fourier coefficients (sine-cosine serie). Due to the different nature of the limit cycles in the states of synchronizationand super-synchronization, we have different modes in the function of interaction in these states.
ACKNOWLEDGMENTS
We gratefully acknowledge financial support by the Spanish Ministerio de Econom´ıa y Competitividad and EuropeanRegional Development Fund under contract RTI2018-097063-B-I00 AEI/FEDER, UE, and by Xunta de Galicia underResearch Grant No. 2018-PG082. Authors are part of the CRETUS Strategic Partnership (AGRUP2015/02). Allthese programs are co-funded by FEDER (UE). [1] S. Strogatz,
Sync: The Emerging Science of Spontaneous Order (Hyperion Press, 2003).[2] J. D. Murray,
Mathematical Biology , 2nd ed. (Springer-Verlag Berlin Heidelberg, 1993).[3] J. Buck, The Quarterly Review of Biology , 265 (1988).[4] C. S. Peskin, Mathematical Aspects of Heart Physiology (Courant Institute of Mathematical Sciences, New York University,New York, 1975).[5] V. Torre, Journal of Theoretical Biology , 55 (1976).[6] R. E. Mirollo and S. H. Strogatz, SIAM Journal on Applied Mathematics , 1645 (1990).[7] A. Ghosh, B. Chance, and E. Pye, Archives of Biochemistry and Biophysics , 319 (1971).[8] E. M. Izhikevich, Dynamical systems in neuroscience: the geometry of excitability and bursting (MIT Press, Cambridge,Mass., 2007). [9] K. Wiesenfeld, P. Colet, and S. H. Strogatz, Physical Review E - Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics , 1563 (1998).[10] A. Hohl, A. Gavrielides, and T. Erneux, Physical Review Letters , 4745 (1997).[11] J. A. Acebr´on, L. L. Bonilla, C. J. P´erez Vicente, F. Ritort, and R. Spigler, Reviews of Modern Physics , 137 (2005).[12] H. Nakao, Contemporary Physics , 188 (2016).[13] I. Z. Kiss, Y. Zhai, and J. L. Hudson, Phys. Rev. Lett. , 248301 (2005).[14] A. F. Taylor, M. R. Tinsley, F. Wang, Z. Huang, and K. Showalter, Science , 614 LP (2009).[15] M. Tinsley, A. Taylor, Z. Huang, F. Wang, and K. Showalter, Physica D: Nonlinear Phenomena , 785 (2010).[16] A. F. Taylor, M. R. Tinsley, F. Wang, and K. Showalter, Angewandte Chemie , 10343 (2011).[17] G. Ghoshal, A. P. Mu˜nuzuri, and J. P´erez-Mercader, Scientific Reports , 19186 (2016).[18] D. A. Prince and B. W. Connors, Advances in neurology , 275 (1986).[19] M. de Curtis and G. Avanzini, Progress in Neurobiology , 541 (2001).[20] T. J. Walker, Science , 891 (1969).[21] B. Pietras and A. Daffertshofer, Physics Reports , 1 (2019).[22] A. Dhooge, W. Govaerts, Y. A. Kuznetsov, H. G. E. Meijer, and B. Sautois,Mathematical and Computer Modelling of Dynamical Systems , 147 (2008).[23] B. Ermentrout, Simulating, Analyzing, and Ani-mating Dynamical Systems: A Guide toXPPAUT for Researchers andStudents (SIAM, Philadelphia, 2002) p. 290.[24] Bo Peng, G. Vilmos, and K. Showalter, Philosophical Transactions of the Royal Society of London. Series A: Physical and Engineering Sciences , 275 (1991).[25] M. Brøns and K. Bar-Eli, Journal of Physical Chemistry , 8706 (1991).[26] M. Krupa and P. Szmolyan, Journal of Differential Equations , 312 (2001).[27] J. Hesse, J.-H. Schleimer, and S. Schreiber, Physical Review E , 52203 (2017).[28] V. Alvarez-Maubecin, F. Garc´ıa-Hern´andez, J. T. Williams, and E. J. Van Bockstaele,The Journal of Neuroscience20