Robust Weak Chimeras in Oscillator Networks with Delayed Linear and Quadratic Interactions

TL;DR

This paper introduces a method to generate robust chimera states in oscillator networks with delayed linear and quadratic interactions, demonstrated in chemical and electrochemical systems, bridging phase models and real-world networks.

Contribution

A general approach using delayed interactions with linear and quadratic terms to produce robust chimeras in oscillator networks, applicable to real-world systems.

Findings

Successfully generated robust chimeras in the Brusselator model.

Validated the approach with electrochemical oscillator experiments.

Bridges gap between phase models and real oscillator networks.

Abstract

We present an approach to generate chimera dynamics (localized frequency synchrony) in oscillator networks with two populations of (at least) two elements using a general method based on delayed interactions with linear and quadratic terms. The coupling design yields robust chimeras through a phase-model-based design of the delay and the ratio of linear and quadratic components of the interactions. We demonstrate the method in the Brusselator model and experiments with electrochemical oscillators. The technique opens the way to directly bridge chimera dynamics in phase models and real-world oscillator networks.

Figures

Equations

$$g(\varphi )=\sin (\varphi -\alpha )+r\sin (2(\varphi -\alpha ))$$

$$\begin{array}{rl}{\dot{\theta }}_{\sigma ,1}& =\omega +g({\theta }_{\sigma ,2}-{\theta }_{\sigma ,1})\\ & +\epsilon \left(g({\theta }_{\kappa ,1}-{\theta }_{\sigma ,1})+g({\theta }_{\kappa ,2}-{\theta }_{\sigma ,1})\right)\end{array}$$

$$\begin{array}{rl}{\dot{\theta }}_{\sigma ,2}& =\omega +g({\theta }_{\sigma ,1}-{\theta }_{\sigma ,2})\\ & +\epsilon \left(g({\theta }_{\kappa ,1}-{\theta }_{\sigma ,2})+g({\theta }_{\kappa ,2}-{\theta }_{\sigma ,2})\right)\end{array}$$

$$\begin{array}{rl}{\dot{\psi }}_1& =g({\psi }_2-{\psi }_1)-g(0)\\ & +\epsilon \left(2g(-{\psi }_1)-g({\psi }_1)-g({\psi }_2)\right)\end{array}$$

$$\begin{array}{rl}{\dot{\psi }}_2& =g({\psi }_1-{\psi }_2)-g(0)\\ & +\epsilon \left(2g(-{\psi }_2)-g({\psi }_1)-g({\psi }_2)\right).\end{array}$$

$${\dot{x}}_{\sigma ,k}$$

$${\dot{y}}_{\sigma ,k}$$

$$g_{\ell }=Z_{-\ell }{h}_{\ell }.$$

$$p_{\sigma ,k}(t)=\sum _{\kappa ,j\in \left\{1,2\right\}}{K}_{\kappa \sigma }h(x_{\kappa ,j}(t-\tau )),$$

$$\begin{array}{rl}h(x_{\kappa ,j}(t))& =k_1\left({\bar{x}}_{\kappa ,j}(t-{\tau }_1)-{\bar{x}}_{\kappa ,j}(t-{\tau }_1-{\tau }_{\text{Br}})\right)\\ & +k_2\left(({\bar{x}}_{\kappa ,j}(t-{\tau }_2{)}^2+{\bar{x}}_{\kappa ,j}(t-{\tau }_2-{\tau }_{\text{Br}}{)}^2\right),\end{array}$$

$$V_{\sigma ,k}(t)=V_0+K\sum _{\kappa ,j\in \left\{1,2\right\}}{K}_{\kappa \sigma }h({\bar{E}}_{\kappa ,j}(t-\tau )),$$

$$\begin{array}{rl}h(E_{\kappa ,j}(t))& =k_1\left({\bar{E}}_{\kappa ,j}(t-{\tau }_1)-{\bar{E}}_{\kappa ,j}(t-{\tau }_1-{\tau }_{\text{Ex}})\right)\\ & +k_2\left(({\bar{E}}_{\kappa ,j}(t-{\tau }_2{)}^2+{\bar{E}}_{\kappa ,j}(t-{\tau }_2-{\tau }_{\text{Ex}}{)}^2\right)\end{array}$$

$$\begin{array}{rl}{\bar{g}}_{\text{Ex}}(\varphi )& =-0.012\cos (\varphi )+0.051\sin (\varphi )\\ & +0.003\cos (2\varphi )-0.021\sin (2\varphi )\end{array}$$

Full text

††preprint: 10.1103/PhysRevLett.119.168301

Phase models provide mathematical descriptions of weakly coupled oscillatory systems. The state of each unit is described by a single variable, its phase, and the effect of coupling is determined by the phase velocity as a function of the phase difference of the coupled elements Kuramoto (1984); Acebrón et al. (2005); Ashwin and Swift (1992). They capture collective dynamical phenomena (e.g., synchronization and dynamical differentiation) of even very large networks of oscillators, as it was demonstrated with electrochemical Kiss et al. (2005) and neural oscillations Hansel et al. (1993); Ashwin et al. (2016), and superconducting Josephson junctions Wiesenfeld et al. (1998). Phase-model-based approaches have also been effective to induce desirable synchronization patterns with external signals, e.g., with desynchronization, and stable and itinerant cluster dynamics Kiss et al. (2007); Rosenblum and Pikovsky (2004); Orosz et al. (2009).

Many biological systems, however, operate at an intermediate level of (frequency) synchronization Uhlhaas and Singer (2006). Collective dynamics where oscillators are only locally frequency synchronized—commonly know as chimeras—are striking examples for the rich dynamics that arise even in identical units Panaggio and Abrams (2015); Schöll (2016) that are of relevance in applications Bick and Martens (2015). Much theoretical effort has focused on chimeras in phase oscillator networks. These range from explicit bifurcation analyses Abrams et al. (2008) to a mathematically rigorous notion of a chimera—a weak chimera is characterized by angular frequency synchronization along trajectories—and corresponding existence results Ashwin and Burylko (2015); Bick and Ashwin (2016); Bick (2017a). At the same time, carefully designed experiments with chimera dynamics have only drawn inspiration from the phase oscillator results Tinsley et al. (2012); Hagerstrom et al. (2012); Martens et al. (2013); Wickramasinghe and Kiss (2013) rather than relate directly to them. Indeed, general experimental realization conditions for robust chimeras (as asymptotic dynamics which arise despite the inherent heterogeneities) are difficult to formulate because of the complexities of the experimental systems. For example, the electrochemical chimera system Wickramasinghe and Kiss (2013) lasted only 100 cycles, required many connections (at least 20 units with 140 connections) and showed chimera dynamics with unrealistically uniform system with natural frequency differences less than $0.1\%$.

In this Letter, we show that very robust chimeras arise in a small oscillator network of only two populations of two elements, when the interactions among the elements are designed in general way with weak linear and quadratic, time-delayed interactions. The interactions are based on a phase model for which we predict the emergence and bifurcations of weak chimeras. The effective design is achieved by the generalization of a feedback approach previously used to induce collective dynamics of globally coupled networks Kiss et al. (2007); Kori et al. (2008); Rusin et al. (2010) to complex network structures. We verify our approach in numerical simulations of the Brusselator model and experiments with electrochemical oscillators to observe weak chimeras in these systems.

*Weak chimeras in networks of phase oscillators—*We consider the dynamics of $M=2$ populations of $N=2$ phase oscillators where the phase interaction between oscillators is determined by the coupling function

[TABLE]

with parameters $\alpha,r\in\mathbb{R}$. More precisely, let the phase $\theta_{\sigma,k}\in\mathbf{T}:=\mathbb{R}/2\pi\mathbb{Z}$ of oscillator $k$ in population $\sigma\in\left\{1,2\right\}$ evolve according to

[TABLE]

where $\kappa=3-\sigma$, $\omega=1$ is the intrinsic frequency of each oscillator 111Note that we can set $\omega=1$ without loss of generality in the phase equations (2) by going into a suitable co-rotating frame. Moreover, in the phase dynamics we omit the global coupling strength $K$ introduced later since we can assume $K=1$ by rescaling time appropriately., and $\varepsilon$ is the interpopulation coupling parameter; see Fig. 1(a) of the Supplemental Material for a sketch of the network topology. If $\theta(t)$ is a trajectory of (2) with initial condition $\theta(0)=\theta^{0}$ then let $\hat{\theta}(t)$ be a continuous lift of $\theta$ to $\mathbb{R}$. With $\Omega_{\sigma,k}(T):=\frac{1}{T}\hat{\theta}_{\sigma,k}(T)$ we have the asymptotic average angular frequency $\Omega_{\sigma,k}=\lim_{T\to\infty}\Omega_{\sigma,k}(T)$ of oscillator $(\sigma,k)$. Recall that the characterizing feature of a weak chimera as a particular invariant set $A\subset\mathbf{T}^{MN}$ is frequency synchrony (and lack thereof): For all trajectories with initial conditions $\theta^{0}\in A$ we have distinct oscillators $(\sigma,k),(\eta,j),(\rho,\ell)$ such that $\Omega_{\sigma,k}=\Omega_{\eta,j}\neq\Omega_{\rho,\ell}$; see Ashwin and Burylko (2015); Bick and Ashwin (2016); Bick (2017a) for a precise definition.

If $\varepsilon=0$, the populations in (2) are uncoupled, which gives rise to invariant subspaces. Each population evolves on $\mathbf{T}^{2}$, and, for a moment, we suppress the population index $\sigma$. The set $\textrm{S}=\left\{\theta_{1}=\theta_{2}\right\}$ corresponds to full phase synchrony and $\textrm{D}=\left\{\theta_{1}=\theta_{2}+\pi\right\}$ denotes the splay phase where oscillators are in antiphase. The asymptotic average frequencies of the oscillators can be written in terms of the coupling function $g$: We have $\Omega_{k}(\theta^{0})=\omega+g(0)$ for $\theta^{0}\in\textrm{S}$ and $\Omega_{k}(\theta^{0})=\omega+g(\pi)$ for $\theta^{0}\in\textrm{D}$. Moreover, $g$ determines the stability of S and D. If $g$ has only a single harmonic, $r=0$, then full synchrony S and D exchange stability at $\alpha=\pm\frac{\pi}{2}$ in a degenerate bifurcation. A second nontrivial harmonic, $r\neq 0$, breaks this degeneracy, that is, for $\alpha\approx\pm\frac{\pi}{2}$ there is a branch of stable (relative) equilibria for $r>0$ Rusin et al. (2010) and a region of bistability between S and D for $r<0$.

For phase shifts $\alpha\approx\frac{\pi}{2}$ and $r<0$, the system (2) now supports weak chimeras for a wide range of parameter values $\varepsilon>0$. Such chimeras arise as perturbations of $\textrm{D}\times\textrm{S}$ 222Since the system is symmetric with respect to permuting the populations there are also weak chimeras close to $\textrm{S}\times\textrm{D}$ for small $\varepsilon$. for small $\varepsilon>0$ Ashwin and Burylko (2015). For the dynamics (2) the space $\mathbf{T}\textrm{S}:=\mathbf{T}^{2}\times\textrm{S}=\left\{\theta_{2,1}=\theta_{2,2}=:\vartheta\right\}$, where the second population is phasesynchronized, is dynamically invariant and its transversal stability is determined by $g^{\prime}(0)$ Ashwin and Swift (1992); Bick and Ashwin (2016). The dynamics on $\mathbf{T}\textrm{S}$ are determined by the phase differences $\psi_{k}:=\theta_{1,k}-\vartheta$ which evolve as

[TABLE]

in a frame that rotates (not necessarily uniformly) with the synchronized population $\theta_{\sigma,k}$. The set $\textrm{S}\textrm{S}:=\textrm{S}\times\textrm{S}=\left\{\psi_{1}=\psi_{2}\right\}\subset\mathbf{T}\textrm{S}$ is dynamically invariant. Fig. 1(a) shows the stable limit cycle that corresponds to the stable weak chimera in the full system and in-phase and antiphase synchronized clusters in the phase plane of (3) for $\alpha=1.57$, $r=-0.3$, and $\varepsilon=0.1$. For initial conditions in $\mathbf{T}\textrm{S}$, trajectories converge either to the weak chimera or to equilibria on SS; thus, we calculated the size of the basin of attraction of the weak chimera periodic orbit as the complement of the basin of attraction $\mathcal{B}(\textrm{S}\textrm{S})$ 333We estimated the size of $\mathcal{B}(\textrm{S}\textrm{S})$ numerically by checking the fraction of whether a trajectory with initial condition on a uniform grid had approached SS after $T=300$ time units up to tolerance $10^{-3}$.. The basin of attraction of the stable periodic orbit shrinks as $\varepsilon\geq 0$ is increased—shown in Fig. 1(b)—and the periodic orbit becomes unstable in a pitchfork bifurcation of limit cycles at $\varepsilon\approx 0.64$. The unstable limit cycle is ultimately destroyed at $\varepsilon\approx 0.715$ in a global bifurcation. The stable weak chimeras here are robust against small perturbations of the system (as hyperbolic limit cycles) and, in contrast to systems where $g$ has only a single harmonic Panaggio et al. (2016), exist for a wide range of parameters with a relatively large basin of attraction. This makes them suitable for realization in limit cycle oscillators through feedback.

*Feedback induces weak chimeras in the Brusselator—*We first illustrate our engineering approach to obtain weak chimeras in limit cycle oscillator systems using the Brusselator model, a simple two-variable ordinary differential equation system that admits a Hopf bifurcation Glansdorff and Prigogine (1971). For real parameters $A$ and $B$, define $f(x,y)=\frac{B}{A}x^{2}+2Axy+x^{2}y$. Let $p_{\sigma,k}(t)$ be a control signal for the $k$th oscillator in population $\sigma\in\left\{1,2\right\}$ whose dynamics are given by

[TABLE]

where $K$ is the total gain for the control signal that is applied to the first component. Fix $A=1$ and $B=2.3$. For $K=0$, each oscillator has a stable limit cycle with angular speed $\omega_{\text{Br}}=0.977$ and period $T_{\text{Br}}=2\tau_{\text{Br}}=2\pi\omega_{\text{Br}}^{-1}$.

For appropriately chosen control and sufficiently small $K$, the network of Brusselator oscillators has a desired phase reduction; see also Kori et al. (2008). More precisely, given a uniformly increasing phase variable $\dot{\phi}_{\sigma,k}=1$ on the limit cycle for $K=0$ and a target interaction function $g(\phi)=\sum_{\ell\in\mathbb{Z}}g_{\ell}\exp(-i\ell\phi)$ for the phase dynamics, then the feedback $h(\phi)=\sum_{\ell\in\mathbb{Z}}h_{\ell}\exp(-i\ell\phi)$ can be obtained from the phase response curve $Z(\phi)=\sum_{\ell\in\mathbb{Z}}Z_{\ell}\exp(-i\ell\phi)$ by solving

[TABLE]

For feedback control, the $h_{\ell}$ may be expressed in terms of the waveform $x(\phi)=\sum_{\ell\in\mathbb{Z}}a_{\ell}\exp(-i\ell\phi)$ which yields a set of equations involving the Fourier coefficients of the waveform, the phase response curve, and the target interaction function $g$.

We realized the network (2)—up to a rescaling of time—with the feedback signal

[TABLE]

where the feedback gains $K_{\kappa\sigma}=1$ if $\kappa=\sigma$ and $K_{\kappa\sigma}=\varepsilon$ otherwise determine the network topology and $\tau$ is a global feedback delay. Let $\bar{x}_{\kappa,j}=x_{\kappa,j}-a_{0}$, where $a_{0}$ is the zeroth Fourier coefficient of $x_{\kappa,j}$ as above (similarly, write $\bar{g}=g-g_{0}$). Now for a given coupling function (1) and second order time-delayed feedback

[TABLE]

the feedback parameters $k_{s}$ and $\tau_{s}$ are readily computed Kiss et al. (2007). Since the delay $\tau_{\text{Br}}$ is equal to half of the oscillations period, the linear feedback term will not have a second harmonic and the quadratic term will not have a first harmonic. Moreover, each delay effectively acts as a phase shift for the coupling function (as long as $K\tau$ is small) Kiss et al. (2007); Kori et al. (2008). Thus, choosing $\tau_{1}$ and $\tau_{2}$ to yield pure $\sin$ interaction for the first- and second-order feedback and then setting $k_{2}/k_{1}=r$ yields the coupling function (1) for $\alpha=0$. For $\alpha\neq 0$, set the global delay to $\tau=\omega_{\text{Br}}^{-1}\alpha$. This strategy yields a good approximation $g_{\text{Br}}$ of the target interaction function (1) as shown in Fig. 2(a) 444The phase response curve was calculated using XPP Ermentrout (2002) to evaluate (5). without the need for extensive nonlinear fitting as in previous approaches Kiss et al. (2007).

Subject to feedback, the network of Brusselator oscillators (4) gives rise to weak chimeras. Fig. 2(b) shows the evolution of $x_{\sigma,k}$ and $\langle\phi_{\sigma,k}\rangle$—the average of $\phi_{\sigma,k}$ over one cycle—for $\varepsilon=0.1$. The average angular frequencies $\Omega_{\sigma,k}$ (calculated from the phase $\phi_{\sigma,k}$ of each oscillator) of the synchronized and antiphase population are distinct. Stable weak chimeras exist for a range of coupling parameters $\varepsilon$ as shown in Fig. 2(c), comparable to the predictions obtained from the phase oscillator dynamics (2). The weak chimeras are robust to adding small variations in the oscillators (not shown).

*Experimental realization in electrochemical oscillators—*We built an experimental system with four (two populations of two) oscillatory chemical reactions that can be coupled through linear and quadratic feedback with a delay. Each oscillatory reaction occurs on the surface of a $1.00\leavevmode\nobreak\ \mathrm{mm}$ diameter nickel wire in $3\leavevmode\nobreak\ \mathrm{M}$ sulfuric acid. Because the disk electrodes are placed far from each other (about $3\leavevmode\nobreak\ \mathrm{mm}$ spacing) and the potential drop in the electrolyte is very small (about $0.1\leavevmode\nobreak\ \mathrm{mV}$), the oscillators do not show synchronization without the presence of additional coupling means Kiss et al. (2002). A multichannel potentiostat interfaced with a real-time Labview controller sets the potential $V_{\sigma,k}(t)$ of the wires individually with respect to a $\textrm{Hg/Hg}_{2}\textrm{SO}_{4}/\textrm{sat.\leavevmode\nobreak\ K}_{2}\textrm{SO}_{4}$ reference electrode. The currents $I_{\sigma,k}(t)$ of the four electrodes and a Pt-coated Ti electrode are recorded and converted to electrode potentials $E_{\sigma,k}(t)=V_{\sigma,k}(t)-I_{\sigma,k}(t)R_{\textrm{ind}}$, where $R_{\textrm{ind}}=1\leavevmode\nobreak\ \mathrm{kOhm}$ is an individual resistance attached to each wire Kiss et al. (2007); see Supplementary Material for more details on the experimental setup. The electrode potential is corrected for offset, $\bar{E}_{\sigma,k}=E_{\sigma,k}-o$, where $o$ is a time averaged electrode potential. For the feedback, the circuit potentials of each wire is set to

[TABLE]

where $K_{\kappa\sigma}$ determines the network topology, $K$ is the total feedback gain, $\tau$ is the global delay, and

[TABLE]

is the feedback as in (7). As above, the delay $\tau_{\text{Ex}}$ is set to be equal to half a period of the uncoupled oscillators: For the experimental setup with no coupling, $K=0$, and potential set to $V_{0}=1160\leavevmode\nobreak\ \mathrm{mV}$ the electrodissolution process is oscillatory with a natural frequency of about $0.45\leavevmode\nobreak\ \mathrm{Hz}$. During the experiment, a 2—3$\leavevmode\nobreak\ \mathrm{mHz}$ difference in natural frequencies between the electrodes was maintained.

Initial trials allowed us to determine the feedback parameters to get the desired coupling functions (1) using the same strategy as in the numerical simulations. Employing pure first- and second-order feedback gains, we set $k_{1}=0.22$, $k_{2}=2.0\leavevmode\nobreak\ \mathrm{V}^{-1}$, $\tau_{1}=\tau_{2}=\tau=0$. With these parameters, we determined the phase interaction function (using a self-feedback method Rusin et al. (2010)). Fig. 3(a) shows that the experimental phase interaction function $g_{\text{Ex}}$ approximates the desired interaction function (1) with $r=-0.4$ and $\alpha=0$ very well. In terms of Fourier coefficients, we obtained

[TABLE]

which shows that there are weak cosine and strong first and second harmonic sinusoidal components with $r=-0.41$. Adding a global delay of $\tau=0.51\leavevmode\nobreak\ \mathrm{s}$, the uncoupled populations, $\varepsilon=0$, exhibited bistability between in-phase and antiphase oscillatory states for $K=0.52$ (see Supplementary Material). This choice of parameters corresponds to a phase shift of $\alpha=1.44$ in the phase model, is expected to fall within the chimera regime and was used in all the following experiments. Before the experiments, the phases of oscillators in populations 1 and 2 were set to anti- and in-phase configurations, respectively, to provide appropriate initial conditions for the weak chimera.

We observed weak chimeras in the experimental setup for a range of coupling parameters $\varepsilon\geq 0$. First note that, if there are no interpopulation connections, $\varepsilon=0$, there is a very large dynamically induced frequency difference of about $18\leavevmode\nobreak\ \mathrm{mHz}$ between population 1 in an antiphase and population 2 in an in-phase configuration; see Figs. 3(c) and 3(d). When the coupling between the populations was increased to $\varepsilon=0.1$, the populations remain approximately in the anti- and in-phase configurations [see Fig. 3(b)] but now exhibited oscillations due to the interaction between populations. Importantly, the two populations exhibited phase drifting behavior relative to each other; this state thus represents a weak chimera state. As is shown in Fig. 3(c), the frequency difference between the populations in the chimera state is much larger (about 9 times) than the frequency difference without interpopulation coupling. We observed a chimera state for a large interpopulation coupling strength up to $\varepsilon=0.5$; see Figs. 3(c) and 3(d). As $\varepsilon$ was increased, the amplitude of the phase difference oscillations of the synchronized population increased. With strong interpopulation coupling at $\varepsilon=0.8$ the weak chimera breaks down and the two populations became phase locked (see Supplementary Material).

*Discussion—*We showed that a simple network of two populations of two elements, coupled through a linear and quadratic amplification with a delay of half the period, can generate very robust chimera patterns with strong phase slipping behavior between the populations. The induced chimeras do not rely on amplitude dynamics (e.g., from chaotic Hart et al. (2016) or amplitude Schmidt et al. (2014) clusters). Similar dynamics are expected with any nonlinear oscillatory system with a phase interaction function that has strong first harmonic and weak second harmonic components. While oscillations very close to a Hopf bifurcation typically have dominant first harmonics in the interaction functions, second harmonics arise naturally away from the Hopf bifurcation point (unless there is some particular symmetry in the individual dynamics) Kori et al. (2014). Therefore, we expect that systems in which the oscillations occur through Hopf bifurcations (e.g., class 2 neurons or resonators Izhikevich (2000)) can exhibit the chimeras. Moreover, interactions between oscillators are often nonlinear in nature; examples include gating mechanisms that limit the transduction of the coupling signals signals, as was demonstrated with glycolytic oscillations Dhumpa et al. (2014) or through synaptic coupling of neurons Kandel et al. (2000). Thus, our results give insights into how chimera dynamics emerge from the first two orders of the interaction, or could be induced if these interactions can be tuned.

Our results directly bridge complex collective dynamics, such as chimeras, in phase models and real-world oscillator networks. On the one hand, this approach elucidates how complex collective dynamics, which arise in phase oscillator networks, translate into networks of limit cycle oscillators. On the other hand, we anticipate insights into the limitations of phase oscillator models to describe collective dynamics of real-world networks beyond clustering; for example, this should further clarify the need to include higher-order interaction terms Rosenblum and Pikovsky (2007); Ashwin and Rodrigues (2016); Bick et al. (2016) in the phase dynamics as they facilitate spatiotemporal dynamics of localized synchrony Bick (2017b).

*Acknowledgements—*C.B. acknowledges the warm hospitality at SLU and helpful discussions with Peter Ashwin. C.B. has received funding from the People Programme (Marie Curie Actions) of the European Union’s Seventh Framework Programme (FP7/2007–2013) under REA Grant Agreement No. 626111. I.Z.K. acknowledges support from National Science Foundation Grant No. CHE-1465013.