Lateral inhibition provides a unifying framework for spatiotemporal pattern formation in media comprising relaxation oscillators

TL;DR

This paper demonstrates that diverse systems of relaxation oscillators with different local kinetics and connection topologies can produce similar spatiotemporal patterns through diffusive coupling via their inactivation component, revealing a universal behavior.

Contribution

It introduces a unifying framework showing how lateral inhibition leads to fundamental pattern classes in relaxation oscillator media, explaining observed universality.

Findings

Identification of two basic pattern classes: synchronized clusters and oscillation death.

Universal emergence of similar patterns across different oscillator systems.

Explanation of pattern formation through lateral inhibition mechanisms.

Abstract

The collective dynamics seen in a wide variety of chemical, biological and ecological systems involve interactions between relaxation oscillators that typically involve fast activation process coupled with a slower inactivation. In this paper, we show that systems of such oscillators having distinct kinetics governing local dynamical behavior and whose interactions are described by different connection topologies, can exhibit strikingly similar spatiotemporal patterns when diffusively coupled via their inactivation component. We explain the apparent universality of this global behavior by showing that relaxation oscillators interacting via lateral inhibition will generally yield two basic classes of patterns, viz., one comprising one or more clusters of synchronized oscillators while the other is a time-invariant spatially inhomogeneous state resulting from oscillation death. All observed collective states can be interpreted either as specific instances of these fundamental patterns or as resulting from their competition.

Figures

Equations

$$\begin{array}{rl}d{u}_i/dt& =\mathcal{F}(u_i,v_i)+D_u{\Sigma }_{j\in S_i}(u_j-u_i)/k_i,\\ d{v}_i/dt& =\mathcal{G}(u_i,v_i)+D_v{\Sigma }_{j\in S_i}(v_j-v_i)/k_i,\end{array}$$

$$\frac{d{x}_i^{(p)}}{dt}=F_{(p)}(x_i^{(1)},x_i^{(2)},\dots ,x_i^{(R)})+D_{(p)}\sum _{j\in S_i}\frac{x_j^{(p)}-x_i^{(p)}}{k_i},$$

$$\begin{array}{rl}\frac{d{u}_i}{dt}& =\mathcal{F}(u_i,v_i)+D_u\sum _{j\in S_i}\frac{u_j-u_i}{k_i},\\ \frac{d{v}_i}{dt}& =\mathcal{G}(u_i,v_i)+D_v\sum _{j\in S_i}\frac{v_j-v_i}{k_i},\end{array}$$

$$\begin{array}{rl}{F}_u(u,v,w)& ={\alpha }_1-{\beta }_{1}u\left(\frac{w^{n_1}}{k_1^{n_1}+w^{n_1}}\right)+{\alpha }_4(1-u)\left(\frac{u^{n_4}}{k_4^{n_4}+u^{n_4}}\right),\\ F_v(u,v,w)& ={\alpha }_2(1-v)\left(\frac{u^{n_2}}{k_2^{n_2}+u^{n_2}}\right)-{\beta }_{2}v,\\ F_w(u,v,w)& ={\alpha }_3(1-w)\left(\frac{v^{n_3}}{k_3^{n_3}+v^{n_3}}\right)-{\beta }_{3}w.\end{array}$$

$$\begin{array}{rl}\mathcal{F}(u,v)& =B+u^{2}v-(1+A)u,\\ \mathcal{G}(u,v)& =Au-u^{2}v,\end{array}$$

$$\begin{array}{rl}\mathcal{F}(u,v)& =ru\left(1-\frac{u}{K}\right)-q\frac{u}{b+u}v,\\ \mathcal{G}(u,v)& =ϵq\frac{u}{b+u}v-dv,\end{array}$$

$$\begin{array}{rl}\mathcal{F}(u,v)& =u(1-u)(u-\alpha )-v,\\ \mathcal{G}(u,v)& =ϵ(ku-v-b),\end{array}$$

Full text

The ubiquity of patterns in nature has, for several decades, stimulated efforts at understanding possible mechanisms that underlie their emergence Cross1993 ; Ball1999 . Such patterns can be manifested in both space and time, with perhaps the most widespread instance of the latter being provided by systems that exhibit relaxation oscillations (see, e.g., Refs. Buzsaki2004 ; Lenz2011 ; PotvinTrottier2016 ; Guzzo2018 ). In the simplest scenarios, these can be understood as resulting from interactions between an activator component and an inhibitory (or inactivation) component that operate at fast and slow time-scales respectively Izhikevich2000 . Such oscillators can, in turn, interact with each other, which can result in non-trivial emergent collective dynamics in systems ranging from the cell to the food web Glass2001 ; Pikovsky2001 . As competition between neighboring elements is a recurring motif in such systems, it is natural to consider the consequences of lateral inhibition Li2015 ; Collier1996 ; Sternberg1988 ; Blakemore1970 in systems of relaxation oscillators.

In this paper we have explored the collective dynamics in a variety of models arising in chemical, biological and ecological contexts. The common thread connecting these diverse systems is that all of them are described by systems of relaxation oscillators coupled to their neighbors through the diffusion of their inactivation components. This may, for instance, arise in spatially extended ecological habitats comprising several patches, with each exhibiting oscillations in predator and prey populations, where only the predator (acting as the inactivation component) can move across neighboring patches, e.g., as in herbivore-vegetation interactions Turchin2003 . Experimental realizations of such systems involving oscillating chemical reactions in microfluidic devices have demonstrated the existence of several striking patterns Toiya2008 . These include anti-phase synchronization as well as spatially heterogeneous time-invariant patterns. While the latter resemble stripes generated by the Turing mechanism Turing1952 , it has been analytically demonstrated using a generic model that these result from spatially patterned oscillation death Singh2013 . Here we show that the patterns seen in these systems can arise in much more general contexts, specificallyinvolving models characterized by different local dynamics and connection topologies, and which describe processes across a wide range of spatial scales. Such “universality” arises from the fact that all patterns generated by such systems can be seen as either specific manifestations of, or arising through interactions between, two basic classes of patterns. Moreover, these fundamental patterns will be observable in any system where lateral inhibition couples oscillators whose dynamics is governed by interactions between components characterized by widely separate time scales.

The dynamical behavior of a spatially extended system of relaxation oscillators that interact over a general connection topology can be described by the following system of ordinary differential equations:

[TABLE]

where $u_{i}$ and $v_{i}$ ($i=1,2,\ldots N$) represent the activation and inactivation components, respectively. The dynamics of an uncoupled node is specified by the functions $\mathcal{F}$ and $\mathcal{G}$, and are governed by parameters whose values are chosen so as to yield oscillations. The diffusion terms represent interactions of each node $i$ with its neighbors (comprising the set $S_{i}$ having size $k_{i}$) through $u$ and $v$ with coupling strengths $D_{u}$ and $D_{v}$ respectively. The net contribution that each node receives through diffusive interactions is normalized by the number of its neighbors ($k_{i}$) so as to make the results comparable across systems with different connection topologies. In this paper we focus on spatiotemporal patterns seen in systems possessing simple regular geometries, viz., one and two-dimensional lattices. Motivated by examples of lateral inhibition in nature, where communication between neighboring regions occurs almost exclusively through the inactivation component (such as the ones mentioned earlier), we consider the case $D_{u}=0$ for all simulations reported here (unless mentioned otherwise). We have explicitly verified that qualitatively similar results are obtained even when this constraint is relaxed (e.g., when $D_{u}=D_{v}$) provided the kinetics of the inactivation component is sufficiently slower than that of the activation component SI . For most results reported here, $N=20$ (for 1-D) and $60^{2}$ (for 2-D), although we have used other values of $N$ to verify that our results are not sensitively dependent on system size. The boundary conditions are taken to be periodic in order to minimize boundary effects. The equations are solved using a variable step stiff solver. Time units in each case are normalized with respect to the time period of an uncoupled oscillator for the corresponding set of parameter values.

We have carried out simulations with different dynamical systems describing the behavior of individual nodes, corresponding to cell-cycle oscillations (adapted from Ferrell et al. Ferrell2011 ), chemical kinetics (viz., the Brusselator model Prigogine1968 ) and predator-prey interactions (belonging to the Rosenzweig-MacArthur class of models Rosenzweig1963 , see Ref. Turchin2003 ). Despite very different expressions for $\mathcal{F}$ and $\mathcal{G}$ representing the intrinsic activity for each of these systems, upon coupling the oscillators exhibit strikingly similar patterns (Fig. 1). Furthermore, the broad nature of the collective dynamics does not appear to depend appreciably on the dimensionality of the lattice, or indeed, its connectivity. The diversity of patterns that can be generated by these models is illustrated in Fig. 1(a) and (b), corresponding to spatiotemporal activity in a ring of cell-cycle oscillators and a chemical system undergoing auto-catalytic reactions that exhibits periodic activity far from equilibrium, respectively. The different patterns shown in (a) correspond to different arrangements of phase relations between neighboring oscillators that are seen for a broad range of system parameters. These include (i) Gradient Synchronization (GS), characterized by a monotonic change in the phase of oscillators along the array (manifested as a propagating front of activity), a special case of which is (ii) Synchronized Oscillations (SO) in which all oscillators are in the same phase, (iii) Anti-Phase Synchronization (APS) where the phase of neighboring oscillators differ by $\pi$, (iv) states corresponding to generalized phase synchronization, (v) Chimera State (CS) comprising co-occurring oscillating and non-oscillating nodes, and (vi) Spatially Patterned Oscillator Death (SPOD) in which the oscillators are arrested at different phases. Rings comprising oscillators described by the Brusselator model shows patterns similar to those seen in (a) as well as more complex ones SI . As shown in (b), these include one or more phase defects, viz., a local phase relation between a set of neighboring oscillators that is distinct from the rest of the lattice, moving across the domain. These relatively exotic patterns occur over a restricted range of the parameter space and are specific to the nature of $\mathcal{F}$ and $\mathcal{G}$.

Even for higher spatial dimensions, we observe that the collective dynamics of systems characterized by different relaxation oscillator models exhibit similar characteristics. Specifically, on comparing the spatiotemporal activity of two-dimensional lattices of coupled Brusselators [Fig. 1(c)] and interacting predator-prey populations described by the Rosenzweig-MacArthur model [Fig. 1(d)], a common set of patterns is observed. These include higher dimensional analogues of GS shown above for 1-D lattices, as well as, generalized APS (where the phase of an oscillator at time $t$ differs by $\pi$ from that of its neighbor at $t+\delta t$, where $\delta t$ is a short time delay). In both of these cases, we observe that the propagating fronts take the form of spiral waves. For stronger diffusive couplings in both systems, we observe CS and SPOD patterns analogous to those reported above for 1-D systems. This suggests that the effective lateral inhibition implemented by diffusion through the inactivation variable is the dominant factor underpinning the patterns that can be seen in these very different systems.

In order to investigate in detail the collective dynamical phenomena common to the systems of relaxation oscillators shown in Fig. 1 we consider a generic model of such oscillators, viz., the Fitzhugh-Nagumo (FHN) model, to describe the dynamics of each node. As in Eqn. (1), each oscillator is described by a fast activation variable $u$ and a slow inactivation variable $v$. Their time-evolution is governed by the functions $\mathcal{F}(u,v)=u(1-u)(u-\alpha)-v$ and $\mathcal{G}(u,v)=\epsilon(ku-v-b)$, where $\alpha=0.139$, $k=0.6$ are parameters specifying the kinetics, $b$ characterizes the asymmetry (related to the ratio of the time the oscillator spends in the high- and low-value branches of the $u$ nullcline) and $\epsilon=10^{-3}$ is the recovery rate. We have verified that the results reported here are not sensitive to small variations in these values. Moreover, introducing different boundary conditions can yield qualitatively similar results Singh2013 .

Fig. 2 (a) shows the three most general patterns that can be observed in a ring of FHN oscillators, viz., GS, generalized APS and SPOD. In addition to these, the system exhibits all other robust patterns that are seen over a wide range of parameter values in the systems described earlier [see Fig. 1 (a-b)] which can be viewed as either special cases (such as SO and APS) or combinations of the three aforementioned general patterns. For instance, CS corresponds to part of the system being in GS while the remainder converges to SPOD.

To quantitatively analyze the robustness of the observed patterns over the ($b,D_{v}$) parameter space, we have estimated the size of their respective basins of attractions from many ($\sim 10^{3}$) realizations [Fig. 2 (b)]. In order to classify the space into distinct pattern regimes we have used the following order parameters. First, SPOD and CS states are distinguished by determining the number of nodes for which the temporal variance of the activation variable, $\sigma^{2}_{t}(u_{i})$ is zero. This allows us to define the number of non-oscillating cells $N_{no}$ that is used to distinguish between SPOD ($N_{no}=N$) and CS ($0<N_{no}<N$). To distinguish between states where all nodes are oscillating, including SO, APS and GS, we obtain the time-average of the variance of the activation variable calculated over the lattice, $\langle\sigma^{2}_{i}(u)\rangle_{t}$, as well as the corresponding quantity for each of the two sublattices comprising alternating sites. The latter are zero for both SO and APS states, which are then distinguished by determining whether $\langle\sigma^{2}_{i}(u)\rangle_{t}$ is zero (SO) or not (APS). If all three time-averaged variances have finite values, the state is characterized as GS if the times at which the activation variable of different nodes reaches the peak value changes monotonically over the lattice, else it is classified as OTHERS. Note that OTHERS includes a large number of diverse complex spatiotemporal patterns including generalized APS. In practice, the different pattern regimes are identified by specifying thresholds on the above order parameters, whose specific values do not affect the qualitative nature of the results. Parameter regions are marked as GS, SO, APS, SPOD, and CS states if they are obtained for $>50\%$ of random initial conditions (i.e., have the largest basin). While changing the system size do not broadly affect the qualitative features, we note that certain patterns such as SO are harder to obtain in larger systems (as the rapid coordination required for exact synchronization is easier to achieve for smaller systems), while the basins of other patterns (such as APS) shrink considerably.

A striking feature of the ($b,D_{v}$) parameter space diagram for the 1-dimensional array of coupled FHN oscillators is that the three general patterns occur at either end of the coupling strength range, with SPOD being found at the higher values of $D_{v}$ while in the lower range we observe GS as well as anti-phase patterns. The anti-phase collective dynamics manifests itself when the individual oscillator limit cycles are highly asymmetric, corresponding to lower values of $b$. Indeed, in the limit of extreme asymmetry where a node remains in one of the branches (slow or fast) for almost the entire duration of its oscillation period, it can be shown analytically that APS is the only stable state for the system Singh2013 . We would also like to note that while SPOD states resemble Turing patterns Turing1952 , the generative mechanism is quite distinct from that of Turing instability and involves the arrest of oscillators into a heterogeneous stationary state, as demonstrated in Ref. Singh2013 . Consistent with the earlier statement that all patterns other than the three general ones (and their special cases) can be seen as combinations thereof, we observe that these patterns (such as CS) occur in the region between the GS and SPOD regimes. Observation of SO at higher values of $D_{v}$ can be interpreted as a result of increased coordination resulting from stronger coupling between the oscillators.

Investigation of the spatiotemporal dynamics in two-dimensional lattices of $N\times N$ coupled FHN oscillators reveals the existence of patterns similar to those seen in Fig. 1 (c-d), including spiral waves, propagating phase defects and CS [Fig. 2 (c)]. The corresponding $(b,D_{v})$ parameter space [Fig. 2 (d)] displays regimes corresponding to SO, APS, CS, SPOD and OTHERS identified using methods similar to those used for 1-dimensional lattices. The strong qualitative similarity with Fig. 2 (b) is visually apparent. Note that SO is seen to occur over a large region of the parameter space to the exclusion of GS as Fig. 2 (d) is for a small lattice (viz., $N=10$). For larger lattices, signals take longer to traverse the domain making global phase coherence less likely, which results in localized phase coordination manifesting as waves (i.e., GS will dominate).

The qualitative similarity of the parameter space diagrams for $1-$ and $2-$dimensional lattices [Fig. 2 (b) and (d)] suggests that the nature of the pattern regimes seen for a system of such coupled oscillators is independent of the dimensionality. To verify this we now consider a mean-field system of globally coupled FHN oscillators which corresponds effectively to the limit of extremely large number of dimensions (Fig. 3). We observe that the parameter space is dominated by essentially two collective dynamical states, viz., Cluster Synchronization (CLSn) comprising in general $n$ oscillator clusters (each cluster being characterized by the common phase of all the constituent nodes) and the temporally invariant Inhomogeneous Steady State (ISS) where the dynamics of each node is arrested to one of two possible values [see inset in Fig. 3]. We would like to point out that all the observed spatiotemporal patterns mentioned earlier can be viewed as instances or combinations of these two fundamental states. In particular, ISS is equivalent to the SPOD state observed in a finite-dimensional lattice. The spatially symmetric SO state where all oscillators have the same phase, and hence belong to a single cluster, corresponds to CLS1. Similarly, the spontaneously broken spatial symmetry APS state comprising two clusters of oscillators that are exactly $\pi$ out of phase, belongs to CLS2. We also observe other CLSn states corresponding to higher values of $n$ in small regions of the parameter space.

Deviating from the mean-field situation by reducing the number of connections per node will result in the emergence of other robust patterns. In particular, it is possible to observe collective states where the number of clusters is equal to the total number of nodes in the system, i.e., CLSN. This will correspond to GS in lattices with finite coordination number. Its occurrence is inversely related to the communication efficiency, i.e., how rapidly signals coordinate activity across the system. This is governed by the diffusion strength, as well as, the number of connections (relative to the system size), and increasing either may result in merging of clusters that could possibly lead to the globally coherent SO (i.e., CLS1) state.

As it is now apparent that the observed patterns in the $1-$ and $2-$dimensional lattices can be understood as instances of the fundamental patterns CLSn, ISS and combinations thereof, we now focus on how these two collective dynamical states compete with each other at the interface of these two regimes in the FHN parameter space. Fig. 4 shows that close to this boundary the system can converge to either one of the two attractors depending on the (randomly chosen) initial condition. This convergence happens after a period of transient activity that resembles diffusion-mediated coarsening phenomena seen, e.g., in binary mixtures Wen1996 ; Mullin2000 . As seen in Fig. 4 (a), a ring of FHN oscillators in this parameter regime exhibit the rapid creation of several domains of varying sizes, each being either in a CLS or SPOD state. Over time some of these domains expand at the expense of others in a process of competitive growth akin to Ostwald ripening Raatke2002 , with the entire system eventually converging to either CLSn or SPOD states [top and bottom panels of Fig. 4 (a), respectively]. The sizes of the basins of attraction for these two states can be discerned from the asymptotic probability density of $f_{osc}$, viz., the fraction of nodes that belong to any of the domains exhibiting oscillations, whose time evolution is shown in Fig. 4 (b). Fig. 4 (c) shows how, as the system approaches the asymptotic state, the number of distinct domains $n^{\rm domain}$ reduces over time through a process of coalescence. Similar coarsening phenomena leading to any of the two fundamental patterns are also seen for two-dimensional lattices of coupled oscillators [Fig. 4 (d)].

To conclude, we have shown that a variety of similar patterns are exhibited by diverse systems having different local dynamics and connection topologies. This can be explained by noting that all of these patterns are either particular manifestations, or arise through interactions between, two fundamental classes of collective dynamical states. These are characterized either by one or more synchronized clusters, or temporally invariant inhomogeneous patterns. As we show, in general, they will arise from the collective dynamics of a system comprising relaxation oscillators that are coupled through the inactivation components, a feature that is common across the systems that we have considered here. While weak interactions typically generate CLSn, stronger coupling yields ISS. This mechanism of pattern formation, distinct from the classic Turing paradigm, is sufficiently generic to have been observed in experimentally realizable settings, such as in coupled electronic circuits Gambuzza2014 . This opens up the possibility of exploring applications for the dynamics reported here, e.g., in the context of computation Menon2014 . Finally, as diffusion is not the only mechanism through which the dynamical components of a spatially extended system interact, it will be of interest to see how the introduction of other processes, such as advection, will affect the nature of the spatiotemporal patterns generated by such systems. This will, for example, be relevant in the context of predator-prey systems embedded in environments subject to hydrodynamic flows Neufeld2001 (e.g., oceanic plankton populations Abraham1998 ; Neufeld2002 ).

I Model description

In the main text we have described the collective dynamics of a system of coupled relaxation oscillators, which are interacting with each other over a specified connection topology (viz., a ring, a two-dimensional torus or a globally coupled scenario). The spatiotemporal activity of such a system of $N$ oscillators, each of which involves $R$ variables, can be described by a model involving $N\times R$ coupled differential equations, viz.,

[TABLE]

where $x^{(p)}_{i}$ ($i=1,2,\ldots,N$, $p=1,2,\ldots,R$) represents the $p$-th component of the $i$-th oscillator. The dynamics of an uncoupled node is specified by the functions $F_{(p)}$ and governed by parameters whose values are chosen so as to yield oscillations. The diffusion terms represent the interactions of each node $i$ with its neighbors (which comprises the set $S_{i}$ of size $k_{i}$) through the different components $x^{(p)}$ with coupling strengths $D_{(p)}$. The net contribution that each node receives through these diffusive interactions is normalized by the number of its neighbors ($k_{i}$) so as to make the results comparable across systems with different connection topologies.

With the exception of one of the models (viz., the cell-cycle model, described below in I.A), the individual oscillators that we have considered are described by only two components that are responsible for activation and inhibition, respectively. These systems can hence be described by the equations:

[TABLE]

which corresponds to Eq. (1) of the main text. In the following sub-sections we describe in detail each of the four models whose collective dynamics is reported in the main text.

I.A Cell-Cycle model

In the first system that we consider, the individual oscillators describe the cell cycle, the periodic sequence of events in a cell resulting in it dividing into two daughter cells, during early embryogenesis in Xenopus laevis, a frog that is native to sub-Saharan Africa. We have adapted the oscillator description used here from a three-component model developed by Ferrell et al. Ferrell2011a which involves interactions between the proteins CDK1, Plk1 and APC. In the course of the cell cycle, CDK1 activates Plk1, which in turn activates the protein APC that subsequently suppresses CDK1.

Representing the concentrations of CDK1, Plk1 and APC by $u$, $v$ and $w$, respectively, the time-evolution of this system is governed by Eqn. 1 with the functions given by:

[TABLE]

Note that this system, with positive and negative feedback loops, behaves like a relaxation oscillator with distinct fast and slow phases. Typical values of the model parameters which generate oscillations are $\alpha_{1}=0.02$, $\alpha_{2}=3$, $\alpha_{3}=3$, $\alpha_{4}=3$, $\beta_{1}=3$, $\beta_{2}=1$, $\beta_{3}=1$, $k_{1}=0.5$, $k_{2}=0.5$, $k_{3}=0.5$, $k_{4}=0.5$, $n_{1}=8$, $n_{2}=8$, $n_{3}=8$ and $n_{4}=8$.

As APC is a relatively large regulatory complex, one may consider interactions through diffusion of this protein to be negligible. We use the fact that Plk1 indirectly suppresses CDK1 by activating APC, to consider it as the inactivation variable in our simulations. We have therefore chosen the diffusion strengths $D_{u}$ and $D_{w}$ to be [math] and considered only diffusive coupling through the variable $v$. The patterns shown in  Fig.1(a) of the main text are obtained by varying the model parameters $\alpha_{1},\alpha_{2},\beta_{1}$ and the coupling strength $D_{v}$.

I.B Brusselator chemical oscillator model

The next system we consider consists of oscillators that describe the far-from-equilibrium behavior of chemical systems in which the concentrations of some reactants exhibit periodic variations. The specific model used is the Brusselator Prigogine1968a , a simplified description of autocatalytic chemical reactions such as the Belousov-Zhabotinsky reaction. The time evolution of the system is described by Eqn. 2 with the functions $\mathcal{F}$ and $\mathcal{G}$ having the following form:

[TABLE]

where $A,B$ are parameters whose values can be appropriately chosen in order to make the system oscillate. For our simulations, we consider $A>A_{\rm hopf}(=1+B^{2})$, i.e., the parameter regime in which individual elements have undergone a Hopf bifurcation to oscillatory behavior. As mentioned in the main text, in experiments done with Belousov-Zhabotinsky reagents in beads that are separate by columns of oil in microfluidic devices, it is known that only the inactivation variable can effectively diffuse through the oil. We have therefore considered $D_{u}=0$. Fig. S1 shows the different patterns that are obtained by varying the values of the parameters $A$, $B$ and $D_{v}$.

I.C Rosenzweig-MacArthur predator-prey model

We next consider a model system whose individual oscillators describe the dynamics arising from interactions between a prey (whose population defines the fast activation variable $u$) and a predator species (whose population is given by the slow inactivation variable $v$). The specific form for the functions are obtained from the Rosenzweig-MacArthur model Rosenzweig1963a ; Turchin2003a , viz.,

[TABLE]

where $r$ is the intrinsic growth rate, $K$ is the carrying capacity of the prey population, $q$ is the maximum predation rate of the predator, $b$ is the half saturation constant and $\epsilon$ and $d$ represent the efficiency and the death rate of the predator, respectively. For our simulations, we have set $r=1$, $K=1$, $q=2$, $d=0.1$ and $\epsilon=0.1$. There are many trophic interactions in which the prey is immobile (e.g., plants) and the predator is able to graze by moving from one patch to another (e.g., herbivore). In such a context, one can set $D_{u}=0$ and vary $D_{v}$.

I.D FitzHugh-Nagumo model

For the bulk of our simulations, we have considered perhaps the most generic model of relaxation oscillators, viz., the Fitzhugh-Nagumo (FHN) model. Here, each oscillator is described by a fast activation variable $u$ and a slow inactivation variable $v$. As discussed in the main text, the time-evolution of this model is governed by the functions:

[TABLE]

where $\alpha=0.139$ and $k=0.6$ are parameters specifying the kinetics, $b$ characterizes the asymmetry (related to the ratio of the time spent in the high- and low-value branches of the $u$ nullcline) and $\epsilon=10^{-3}$ is the recovery rate. Fig. S2 (a) shows that in the absence of diffusive coupling with its neighbors, oscillation is seen over a range of values of the parameter $b$. Remaining within this interval, we have varied both $b$ and the coupling strength $D_{v}$ to observe a range of collective dynamics [Fig. S2 (b)].

II Coarsening

In the main text we have shown that competition between attractors at the boundary of the GS and SPOD regimes for a ring of FHN oscillators results in coarsening-like dynamics. Beginning with many small domains, each consisting of exclusively oscillating or non-oscillating elements, the collective dynamics exhibits a reduction in the number of domains over time through a process of coalescence. Here we show how the coarsening depends on the location of the system in the $D_{v}-b$ parameter space by gradually increasing the coupling strength $D_{v}$ (keeping the value of $b$ fixed at $0.175$). In Fig. S3, starting from deep within the GS regime [panel (a)] and terminating deep in the SPOD regime [panel(e)], we show spatiotemporal evolution of the inactivation variable $v$. We see that, for lower values of diffusion strength $D_{v}$, gradient synchronization (GS) occurs almost immediately with negligible transients and hence, no coarsening is seen for this regime [panel (a)]. However, on increasing $D_{v}$ coarsening phenomena is observed in the transient period leading to the asymptotic state (either GS or SPOD). Indeed, we note that coarsening is seen over a range of $D_{v}$, with a lower bound is $\sim 1.76\times 10^{-3}$ and upper bound is $\sim 2\times 10^{-3}$. Panels (b-d) show coarsening to GS, extremely long-lived transients and SPOD for the same value of $D_{v}$ that falls in the above range. Panel (e) shows that for high $D_{v}$ convergence to SPOD occurs rapidly (without an appreciable period of transients) for high $D_{v}$.

III Robustness of Collective Dynamics: Results for $\bf{D_{u}=D_{v}}$

For the simulations reported in the main text, we have considered $D_{u}=0$. However, as shown in Fig. S4, we obtain qualitatively similar results even for the case $D_{u}=D_{v}$, i.e., when both activation and inactivation variables are diffusively coupled with the same strength to their neighboring sites.