Dynamical order, disorder and propagating defects in homogeneous system of relaxation oscillators
Reaction-diffusion (RD) mechanisms in chemical and biological systems can yield a variety of patterns that may be functionally important. We show that diffusive coupling through the inactivating component in a generic model of coupled relaxation oscillators give rise to a wide range of spatio-temporal phenomena. Apart from analytically explaining the genesis of anti-phase synchronization and spatially patterned oscillatory death regimes in the model system, we report the existence of a chimera state, characterized by spatial co-occurrence of patches with distinct dynamics. We also observe propagating phase defects in both one- and two-dimensional media resembling persistent structures in cellular automata, whose interactions may be used for computation in RD systems.
Understanding how patterns develop in chemical and biological systems, e.g., in the course of morphogenesis Morelli12 (), is one of the enduring challenges of modern science Cross1993 (). Investigating pattern formation in systems of coupled biochemical oscillators is a promising approach towards this objective and has been applied to study, for instance, the temporal organization of gene activity in cells during the course of development Goodwin64 (). The patterns seen in such systems are the dynamical attractors that can change with the underlying system parameters as the environment evolves with time. Such phenomena have been sought to be explained by reaction-diffusion models that exhibit patterns such as stripes and spots under specific conditions, viz., when an inhibitory chemical species diffuses faster than an activating species, by destabilizing the homogeneous steady state Kondo10 (). Indeed, generalizations of such processes involving differential excitatory and inhibitory interactions between elements is seen in a wide variety of complex systems Meinhardt72 (); Gong12 (); Butler12 (). Recently such mechanisms have also been proposed as a possible basis for computation in biological and chemical systems Liang1995 (); Epstein2007 ().
Time-invariant patterns that are seen in the models mentioned above are, however, only a small segment of the variety seen in nature, many of which exhibit periodic activity. Thus, extending the concept of reaction-diffusion mechanism to systems of interacting relaxation oscillators may allow investigation of spatio-temporal patterns in biological systems, where oscillations are observed across many spatial and temporal scales, e.g., from the periodic variations of intracellular molecular concentrations Orchard83 () to changes in the activity levels of different brain areas Buzsaki04 (). Coherent dynamics of such oscillators can produce functionally important collective behavior such as synchronization Kiss02 () giving rise to different biological rhythms Singh12 (). However, synchronized oscillations is only one of a number of possible collective phenomena that can emerge from such interactions. For example, a recent set of experiments on coupled chemical oscillators in a microfluidic device Toiya08 () have shown that anti-phase synchronization as well as spatially heterogeneous oscillator death states BarEli85 () can occur under different conditions. Extending the mechanism of coupling by lateral inhibition (e.g., via a rapidly diffusing inhibitory chemical species) to arrays of coupled relaxation oscillators, used for modeling biological periodic activity, may reveal the underlying mechanism for a variety of spatiotemporal phenomena seen in real systems.
In this paper, we have analyzed in detail a model of generic relaxation oscillators (each comprising activator and inactivating components) coupled to their nearest neighbors through lateral inhibition via diffusion of the inactivating component. The model reproduces all the spatiotemporal patterns observed in the experiments mentioned earlier and its simplicity allows an analytical understanding of their genesis. In particular, we have given possibly the simplest theoretical demonstration of the existence and stability of an anti-phase synchronized state for coupled relaxation oscillators. In addition to the patterns reported in experiments, we also observe other states such as, attractors corresponding to spatially co-existing dynamically distinct configurations which we term chimera states. Although homogeneous arrays of generic relaxation oscillators have been studied extensively, the observation of such spatially heterogeneous attractors for these systems is a novel finding of our paper. We have performed for the first time a systematic characterization of the basins of attraction for various patterns seen in the model that also demonstrates the unexpected robustness of the chimera states and suggest that all the observed states can be obtained in suitable experiments. In addition, we report the occurrence of phase defect-like discontinuities moving ballistically through the system that on collision with each other can produce complex patterns. We observe analogous structures in two-dimensional media that have a striking resemblance to persistent configurations in cellular automata (CA), e.g., “gliders” in the “Game of Life” CA Conway (), which have been linked to the universal computation capabilities of such systems Wolfram (); Cook04 (). The observation of these patterns in our model system is remarkable considering the simplicity of the underlying dynamics.
The system we have analyzed comprises relaxation oscillators interacting with each other in a specific topology. For the dynamics of individual relaxation oscillators we have used the phenomenological FitzHugh-Nagumo (FHN) equations, which is a generic model for such systems. Each oscillator is described by a fast activation variable and a slow inactivation variable :
where, , are parameters describing the kinetics, characterizes the recovery rate of the medium and is a measure of the asymmetry of the oscillator (measured by the ratio of the time spent by the oscillator at high and low value branches of ). Parameter values have been chosen such that the system is in oscillatory regime. Small variations in the values do not affect our results qualitatively. To investigate spatial patterns generated by interaction between the oscillators, they are arranged in a 1-dimensional chain [Fig. 1 (a)]. In actual chemical experiments, the beads containing the reactive solution are suspended in a chemically inert medium which allows passage of only the inhibitory chemical species Toiya08 (). In our model, the oscillators are diffusively coupled via the inactivation variable . The boundary conditions for the chain take into account the inert medium by including non-reactive passive elements at each end that are diffusively coupled to the neighboring oscillators. The inert medium between the oscillators are not considered explicitly, their volume being relatively small compared to the reservoirs at the boundary. We have verified that inclusion of intermediate non-reactive cells diffusively coupling each pair of oscillators do not affect the fixed-point equilibria of the system or their stability, subject to suitable rescaling of the diffusion constant. The dynamics of the resulting system is described by
where and the diffusion constant represents the strength of coupling between neighboring relaxation oscillators through their inactivation variables. For most results reported here although we have used higher values of upto to verify that the qualitative results are not sensitively dependent on system size. We have also verified that the boundary conditions do not sensitively affect the results by considering periodic boundaries and observing patterns qualitatively identical to those reported here. The dynamical equations are solved using an adaptive Runge-Kutta scheme (e.g., rkf45 GSL ()). The behavior of the system for each set of parameter values and is analyzed over many () initial conditions, with each oscillator having a random phase chosen from a uniform distribution.
Fig. 1 (b-e) shows a variety of asymptotic spatio-temporal patterns that are observed in the model system: synchronized oscillations (SO) with all elements (except those at the boundary) having the same phase (b), anti-phase synchronization (APS) with neighboring elements in opposite phase (c), Spatially Patterned Oscillator Death (SPOD) regime where the oscillators are arrested in various time-invariant states (d) and Chimera States (CS) where oscillating regions co-exist with patches showing negligible temporal variation (e). However, these do not exhaust the range of possible spatio-temporal phenomena that are observed including propagating structures that are discussed later. Note that, both APS and SPOD states have been observed experimentally in chemical systems. Although the latter have sometimes been referred to as “Turing patterns” in the literature, SPOD is distinct as it is not obtained through destabilization of a homogeneous equilibrium (Turing instability) but through a process of oscillator death.
To investigate the robustness of the observed patterns in detail, we numerically estimate the size of their basins of attraction in the () parameter space (Fig. 2). This provides an indication as to whether a pattern can be observed in actual experiments as for this they must be immune to the unavoidable noise present in any experiment. In order to segregate the distinct pattern regimes in the () space [Fig. 2 (a)] we introduce the following order parameters. The number of non-oscillating cells in the bulk of the system, , i.e., cells for which the variance of the activation variable , , is zero, is used to characterize the SPOD () and CS regimes (). The SO and APS states both have all elements in the bulk oscillating. However, SO is distinguished by having all oscillators in the same phase as measured by the variance of the activation variables , = 0, where represents time average. We can also define the synchronization among the oscillators in two distinct sub-lattices, namely, those which occur at even-numbered and those which occur at odd-numbered sites of the chain, as measured by the time-averaged variance of the activation variable, viz., and . This pair of order parameters are zero for both SO and APS states; however, if , it signifies the APS regime. In practice, the different regimes are characterized by thresholds whose specific values do not affect the qualitative nature of the results. Fig. 2 (a) indicates the parameter regions where SO, APS, SPOD and CS states are observed for more than of initial conditions (i.e., they have the largest basin). As already mentioned, the system also exhibits other regimes apart from the above ones, which occur in regions of () parameter space shown in white.
While diffusive coupling in a homogeneous system of oscillators is expected to promote the SO state Kurths (), a striking result from this phase diagram is that for certain parameter values the APS state has a very large basin of attraction [Fig. 2(b)]. The existence of APS is somewhat counter-intuitive as for diffusively coupled identical isochronous oscillators the only stable attractors are synchronized oscillations or oscillator death Kurths (). While anti-phase synchronization has been seen earlier in a pair of identical oscillators Rand80 (), it is not obvious that APS will have a large basin of attraction for an array of oscillators. To understand the origin of such anti-phase oscillations we consider a simple model that captures the essence of relaxation oscillation phenomena and can be solved exactly.
We consider the relaxation limit ( in FHN system) and extreme asymmetry where the limit cycle has a slow segment in which the system spends the entire duration of the oscillation period (the remaining segment of the cycle being traversed extremely fast). Under these considerations, we obtain the one-dimensional dynamical system: where parameterizes the slow part of the limit cycle and can be redefined to belong to the interval . Fig. 3 (a) shows a schematic diagram of the trajectory of the limit cycle, where the system spends almost its entire oscillation period on the solid branch (the return from to , shown by the broken line, is considered to be instantaneous). The model can be exactly solved if is a constant (, say), although the geometrical argument is valid for any arbitrary positive definite function defined over the interval . By appropriate choice of time scale, we can set the period without loss of generality. A system of two such diffusively coupled oscillators can be described by
Given the values of at some arbitrary initial time , the solution of Eqn. (3) at a later time follow the relations:
till time when reaches . After this the larger of is mapped back to (because of the instantaneous nature of the remaining segment of the limit cycle) and in Eqn. (4) is replaced by . Using the above exact solution of the coupled system (3), its Poincare map is constructed in two steps. First, if starts at and starts at some point , we find the location of at some time when (which is then immediately mapped to ). Now, starting with and , when we find the location of : . Using Eqn. (4), with , , and , we get
where is the Lambert W-function. Fig. 3 (b) shows the Poincare map for different values of the coupling strength . The map has one stable and one unstable fixed point, which correspond to the anti-phase synchronized (APS) and synchronized oscillating (SO) states, respectively. Thus, for the model (3) we find that APS state not only exists but is the only stable state. Relaxing the extremal conditions under which this was derived may allow a stable SO state to coexist with the stable APS state Izhikevich_2000 (). The derivation we have shown here is possibly the simplest one exposing the fundamental mechanism for generating APS states in any system of diffusively coupled oscillators that can exhibit anti-phase oscillations.
When the coupling between oscillators in the array is increased to very high values, we observe that the oscillatory regimes (e.g., SO and APS) are replaced by time-invariant spatial patterns such as SPOD (Fig. 2). To understand the genesis of SPOD at strong coupling, we can again focus on a pair of coupled relaxation oscillators in the relaxation limit (). The parameter is chosen such that the -nullcline is placed symmetrically between the two branches of the -nullcline with the oscillator taking equal time to traverse each branch [Fig. 3 (d)]. When the two oscillators (1 and 2) are in opposite branches (as shown in the schematic diagram), the two opposing forces acting on each oscillator, corresponding to the coupling  and the intrinsic kinetics () respectively, can exactly cancel when the coupling is strong resulting in oscillator death. Note that only the component of parallel to the intrinsic force () needs to balance as has no effect in the relaxation limit. The symmetry of the oscillator ensures that the force due to the intrinsic kinetics () for the two oscillators are identical but oppositely directed in the steady state. The occurrence and stabilization of this heterogeneous time-invariant state (as 1 and 2 need to be in different branches, corresponding to very different values of ) is the key to the occurrence of SPOD at strong coupling. At intermediate values of coupling between oscillators in a large array, the competition of this mechanism with the natural oscillatory dynamics (seen at low coupling) may give rise to chimera states that exhibit characteristics of both, i.e., containing patches of elements that are oscillating as well as segments that show SPOD-like pattern. This CS regime is especially interesting as the system exhibits a strikingly heterogeneous dynamical state in spite of the bulk of the array being homogeneous with identical elements coupled in the same fashion. We have explicitly verified that the occurrence of CS is not dependent on boundary conditions by reproducing it also in systems with periodic boundaries. The observation of such states in a generic model of relaxation oscillators suggests that they should be present in a wide class of systems, and indeed similar phenomena have been recently reported in a specific chemical system model Vanag11 (). Note that, the chimera state described here comprises regions with dynamically distinct behavior, and is different from its namesake that refers to co-occurrence of coherent and non-coherent domains Strogatz2004 ().
Apart from the spatio-temporal patterns shown in Fig. 1 (b-e) we also observe attractors corresponding to point-like “phase defects” (i.e., there is a discontinuity of phase along the oscillator array at this point) moving in the background of system-wide oscillations. As seen from a typical example of such states [Fig. 4 (a)], after initial transients these defects move in the medium with interactions between two such entities resulting in either the two getting deflected in opposite directions, or either one or both getting annihilated (unlike defects in non-oscillatory media such as domain walls which annihilate on collision ChaikinLubensky ()). While the boundary for systems with passive elements at the ends is a source of new defects entering the medium, similar persistent structures are also seen in systems with periodic boundary conditions where, in the steady state, a conserved number of defects can reflect off each other indefinitely [Fig. 4 (b)].
To observe how these propagating defects manifest in higher dimensional systems, we consider a 2-dimensional array of coupled oscillators defined on a torus [Fig. 4 (c-d)]. The system can have extremely complicated transient phenomena, and for simplicity we show here only its asymptotic behavior. For a square lattice, we observe that there is a specific configuration of four sites that persists indefinitely (reminiscent of the glider configurations in the 2-dimensional CA “Game of Life” Conway ()) and can move in horizontal or vertical directions. Interaction of such “gliders” can produce complex spatio-temporal patterns. Fig. 4 (d) shows two “gliders” that on collision move off in directions perpendicular to their incident ones supp ().
To conclude, we have shown that a simple model of relaxation oscillators coupled via lateral inhibition can exhibit a wide variety of striking spatio-temporal patterns. Our comprehensive investigation has revealed the global features of the dynamics and robustness of the attractors. As our results are based on a very generic model, it suggests that the patterns may be observed in a range of experimental realizations, such as electronic circuits implementing relaxation oscillators Scholl2010 () and Pt wire undergoing CO oxidation where the system is in an oscillatory regime Ertl (), apart from microfluidic chemical systems mentioned earlier. Our initial exploration of propagating configurations in 2-dimensional media suggests that systems of higher dimensions may yield even more striking discoveries. It is intriguing to explore the possibility of using the propagating defects for computation Liang1995 (); Epstein2007 (), as analogous entities have been used to construct logic gates in CA Conway (). This may tie together two of the groundbreaking ideas of Alan Turing, viz., pattern formation through reaction-diffusion mechanism and universal computation Turing ().
We thank Bulbul Chakraborty and Sudeshna Sinha for helpful discussions and the HPC facility at IMSc for providing computational resources.
- (1) L. G. Morelli, K. Uriu, S. Ares, A. C. Oates, Science 336, 187 (2012).
- (2) M. C. Cross and P. C. Hohenberg, Rev. Mod. Phys. 65, 851 (1993).
- (3) B. C. Goodwin, Symp. Soc. Exp. Biol. 18, 301 (1964).
- (4) S. Kondo and T. Miura, Science 329, 1616 (2010).
- (5) A. Gierer and H. Meinhardt, Kybernetik 12, 30 (1972).
- (6) Z. Gong et al, Proc. Natl. Acad. Sci. USA 109, E234 (2012);
- (7) T. Butleret al, Proc. Natl. Acad. Sci. USA 109, 606 (2012).
- (8) P. Liang, Phys. Rev. Lett. 75, 1863 (1995).
- (9) I.R. Epstein, Science 315, 775 (2007).
- (10) C. H. Orchard, D. A. Eisner and D. G. Allen, Nature 304, 735 (1983); A. Goldbeter, Biochemical Oscillations and Cellular Rhythms (Cambridge Univ. Press, Cambridge, 1997).
- (11) G. Buzsáki and A. Draguhn, Science 304, 1926 (2004).
- (12) I. Z. Kiss, Y. Zhai and J. L. Hudson, Science 296, 1676 (2002).
- (13) R. Singh et al, Phys. Rev. Lett. 108, 068102 (2012).
- (14) M. Toiya, V.K. Vanag, and I.R. Epstein, Angew. Chem. 47, 40, 7753-5, (2008); J. Delgado et al, Soft Matter 7, 3155 (2011).
- (15) K. Bar-Eli, Physica D 14, 242 (1985).
- (16) E.R. Berlekamp, J.H. Conway and R.K. Guy, Winning Ways for your Mathematical Plays (Academic, New York, 1982), Vol. 2.
- (17) S. Wolfram, Rev. Mod. Phys. 55, 601 (1983); S. Wolfram, Phys. Rev. Lett. 54, 735 (1985).
- (18) M. Cook, Complex Systems 15, 1 (2004).
- (19) Runge-Kutta-Fehlberg (4,5) method, GNU Scientific Library (http://www.gnu.org/software/gsl/).
- (20) A. Pikovsky, M. Rosenblum and J. Kurths, Synchronization (Cambridge University Press, Cambridge, 2003).
- (21) R. Rand and P. Holmes, Int. J. Non. Lin. Mech. 15, 387 (1980).
- (22) E.M. Izhikevich, SIAM J. App. Math. 60, 1789 (2000).
- (23) V. K. Vanag, and I. R. Epstein, Phys. Rev. E 84, 066209 (2011).
- (24) D.M. Abrams and S.H. Strogatz, Phys. Rev. Lett. 93, 17, 174102 (2004).
- (25) P.M. Chaikin, T.C. Lubensky, Principles of Condensed Matter Physics (Cambridge Univ. Press, Cambridge, 2000).
- (26) See supplementary material for movies.
- (27) M. Heinrich et al, New J. Phys. 12, 113030 (2010).
- (28) G. Ertl, Science 254, 1750 (1991).
- (29) A. M. Turing, Phil. Trans. R. Soc. Lond. B 237, 641, 37 (1952); A. M. Turing, Mathematical Logic (Eds. R. O. Gandy and C. E. M. Yates, Elsevier, Amsterdam, 2001).