Dynamics of Reaction-Diffusion Oscillators in Star Networks

TL;DR

This study explores the complex dynamics of inhibitory reaction-diffusion oscillators on star networks, revealing three attractors and analyzing how network structure influences phase synchronization and oscillation behavior.

Contribution

It provides experimental and theoretical insights into star network oscillator dynamics, identifying three attractors and quantifying phase shifts related to network topology.

Findings

Identified three dynamical attractors: unlocked, locked, and center silent.

Discovered a logarithmic relationship between phase shift and star degree.

Differentiated effects of network topology and chemistry on phase shifts.

Abstract

We experimentally and theoretically investigate the dynamics of inhibitory coupled self-driven oscillators on a star network in which a single central hub node is connected to $k$ peripheral arm nodes. The system consists of water-in-oil Belousov-Zhabotinsky $\sim$100$\mu$m emulsion drops contained in storage wells etched in silicon wafers. We observed three dynamical attractors by varying the number of arms in the star graph and the coupling strength; ($i$) unlocked; uncorrelated phase shifts between all oscillators, ($ii$) locked; arm-hubs synchronized in-phase with a $k$-dependent phase shift between the arm and central hub, and ($iii$) center silent; central hub stopped oscillating and the arm-hubs oscillated without synchrony. We compare experiment to theory. For case ($ii$), we identified a logarithmic dependence of the phase shift on star degree, and were able to discriminate between contributions to the phase shift arising from star topology and oscillator chemistry.

Figures

Equations

$${\dot{c}}_i^m=F_i^m\left({\mathbf{c}}_i\right)+\sum _{j=1}^N{\mu }_{ij}^m{\mathcal{A}}_{ij}\left(c_j^m-c_i^m\right),$$

$${\dot{\Phi }}_{0a}=k{\mathcal{H}}_{0a}\left(-{\Phi }_{0a}\right)-{\mathcal{H}}_{a0}\left({\Phi }_{0a}\right)=\mathcal{G}\left({\Phi }_{0a}\right).$$

$${\Phi }_{0a}^{\ast }=\pi +\log \left(k\right)\frac{1}{2}{\mathcal{H}}_{0a}\left({\Phi }_{a0}\right){{\left(\frac{d{\mathcal{H}}_{0a}}{d{\Phi }_{a0}}\right)}^{-1}|}_{{\Phi }_{a0}=\pi },$$

Full text

Discrete networks of self-driven oscillators represent a broad class of physical systems Pecora et al. (2014); Sorrentino et al. (2016); Kopell and Ermentrout (1988); Edwards (2005); Reimann et al. (2017); Néda et al. (2000); Pantaleone (2002); Dörfler et al. (2013); Pikovsky et al. (2003); Strogatz (2000); Kuramoto (2012). Idealized models have the virtue of being simple enough to explain network dynamics, but raise the question of whether such models can quantitatively characterize natural systems. Addressing this question has motivated the development of controlled, model experimental systems that operate with varying degrees of autonomy Jia and Kiss (2012); Wickramasinghe and Kiss (2013, 2014); Horvath et al. (2015); Li et al. (2014, 2015); Tompkins et al. (2015); Litschel et al. (2018); Blaha et al. (2019). A particularly important goal of oscillator network research is to identify how (i) individual oscillator’s dynamics, (ii) network topology and (iii) coupling type conspire to create emergent spatio-temporal patterns.

To address this goal, we experimentally study and theoretically model a minimally complex microfluidically-assembled discrete network of Belousov-Zhabotinsky (BZ) chemical oscillators. Notably, our BZ system uses natural physical-chemical processes to produce the oscillator dynamics, network topology and coupling, in contrast to systems in which the oscillators, coupling, or both are mediated through electronic hardware Miyazaki and Kinoshita (2006); Kiss et al. (2006); Pecora et al. (2014); Sorrentino et al. (2016); Totz et al. (2018); Blaha et al. (2019). Thus, our system also demonstrates the potential of engineering stand-alone soft materials to exhibit self-organized spatiotemporal patterns. An application of such materials is to make autonomous, soft robots that run purely on chemicals like living organisms, rather than powered by motors and controlled by computers. In this scenario, the BZ network would function as a central pattern generator found in the nervous system of many animals Litschel et al. (2018). The BZ reaction requires only a few reagents and networks can be fabricated and parallelized through microfluidic techniques Tompkins et al. (2015). The theoretical foundation of BZ chemistry is strong, enabling near-quantitative modeling of both isolated and coupled oscillators. We examine inhibitory coupling, which promotes symmetry-breaking phenomena by preferring anti-phase synchrony between neighbors Vanag and Epstein (2009, 2011); Tompkins et al. (2014); Li et al. (2014); Guzowski et al. (2015); Torbensen et al. (2017); Litschel et al. (2018). We arrange the cells in a network with star topology, consisting of a central hub node connected by $k$-arms to other nodes, as illustrated in Fig. 1. Star networks have been considered in idealized theoretical and experimental studies Pecora and Carroll (1998); Bergner et al. (2012) and are a naturally occurring motif in neural networks that perform cognitive Bonifazi et al. (2009); van den Heuvel and Sporns (2013); Bertolero et al. (2018) and sensorial functions Khaledi-Nasab et al. (2018).

We probe the dynamical states of the system as a function of inter-nodal coupling strength and star degree, thereby focusing on the impact of physical changes to the star network while taking the BZ chemistry as a known in our experiments and models. We compare observations to two levels of theory: a discrete reaction-diffusion model with Vanag-Epstein Vanag and Epstein (2009) chemical dynamics taking place at each node, and a phase model constructed from it that we use to examine the topology-dependence of locking angles. Strikingly, we find that the locking angle between arm and hub nodes has a logarithmic dependence on the number of arm nodes in the limit of weak coupling. The prefactor for this dependence depends only on the interaction function, which can be numerically derived for any oscillator from a model or an experimentally acquired phase response curve Stankovski et al. (2017). Our result, general to all star networks, disentangles oscillator and coupling physics from topological effects by compactly showing how each come together to produce topology-dependent phase locking.

Experimental System. Surfactant stabilized emulsions of 100 micron diameter drops containing the aqueous BZ solution were generated in a fluorinated oil Delgado et al. (2011). A drop of the concentrated emulsion was pipetted onto the etched silicon wafer, a cover glass was laid on top and clamped together, thereby squeezing the BZ drops into the etched network and sealing the device. Each wafer contained hundreds of star networks. The design concept was to make a device that is simple and quick to load, intolerant to failure, and reusable. Typically 90% of the networks fail to load correctly. However, 10 - 20 successes are enough to accumulate statistics. Both the silicon and glass are completely impermeable to all chemical species ensuring that each arm-drop only communicates with the hub-drop Tompkins et al. (2015). The cavities containing the BZ drops are connected with channels designed to be too narrow to house drops, but contain oil, and therefore function as diffusive conduits, as illustrated in Fig. 1B. The BZ reaction oscillates between a reduced and oxidized state of the catalyst. The duration of the oxidized, or activated state, is brief and during this interval a large amount of the inhibitor, bromine, is also generated. The inhibitory coupling between two drops is provided by bromine, which due to its low polarizability readily partitions into the oil Li et al. (2014). All the chemicals, their concentrations, and conditions for producing the emulsion are described in the supplement.

Theory and Model. In a reaction-diffusion network consisting of $i=1\cdots N$ nodes and $m=1\cdots M$ species, the dynamics of the concentration $c_{i}^{m}$ are governed by,

[TABLE]

where $F_{i}^{m}\left(\mathbf{c}_{i}\right)$ models intra-nodal reactions given by the Vanag-Epstein model of BZ chemistry Vanag and Epstein (2009) and the second term captures inter-nodal diffusive transport proportional to the species-dependent coefficient $\mu_{ij}^{m}$; full expressions are found in the supplement BZS . This model assumes that all the chemistry occurs at discrete points, i.e. it ignores concentration gradients within a drop. The point approximation is justified because the width of the oxidation front in the BZ reaction is larger than the size of the droplets, thus each reactor oxidizes uniformly.

Modeling the coupling as linearly proportional to the concentration difference between connected nodes is equivalent to ignoring any chemical reactions and accumulation of chemicals in the oil separating drops, which is justified when the inter-drop gap is much less than the diffusion length scale $l\sim\sqrt{DT}$, where $D$ is the diffusivity of bromine in oil and the timescale is given by the oscillation period $T$ Li et al. (2014, 2015); Tompkins et al. (2015). We consider gaps no greater than 60 $\mu$m. With diffusivity $D\sim\mathscr{O}\left(10^{-9}\right)$ m2s*-1* and oscillation period $T\sim\mathscr{O}\left(10^{2}\right)$s, the diffusion length scale is $l\sim$300$\mu$m; thus we are safely in the quasi-steady regime. Network connectivity is given by the adjacency matrix $\mathbf{\mathcal{A}}$ such that $\mathcal{A}_{ij}=1$ if $i$ and $j$ are connected and $\mathcal{A}_{ij}=0$ otherwise. We restrict our attention to the dynamics of Eqn. 1 on star graphs with $k$ arm nodes, illustrated in Fig. 1.

We additionally employ the method of phase-reduction to create a further simplified model of the phase-locked dynamics BZS ; Izhikevich (2015); Monga et al. (2019). The approach assumes that the state of an isolated multi-variable chemical oscillator with limit-cycle dynamics can be fully described with a single phase variable. For weakly-coupled oscillators, the slow dynamics of the network, obtained by averaging over one oscillation period, take the form $\dot{\Phi}_{i}=\omega+\sum_{j=1}^{N}{\mathcal{A}_{ij}\mathscr{H}_{ij}}\left(\Phi_{ji}\right)$ where $\Phi_{i}$ and $\dot{\Phi}_{i}$ are the phase and instantaneous frequency of the $i^{\mathrm{th}}$ oscillator, respectively, $\omega$ is the frequency of an isolated oscillator (assumed to be identical for all oscillators), $\Phi_{ji}=\Phi_{j}-\Phi_{i}$ is the phase difference between two oscillators, and $\mathscr{H}_{ij}$ is a numerically calculated interaction function describing how the instantaneous frequency of oscillator $i$ is changed by the presence of oscillator $j$, derived from Eqn. 1 and shown in Fig. 2 111We note that by conventionIzhikevich (2015) $\mathscr{H}_{ij}$ is defined to be a function of $\Phi_{j}-\Phi_{i}$; however, in Fig. 2 we plot as a function of $\Phi_{i}-\Phi_{j}$ to compare to $\mathscr{G}$. This change in variables mirrors the function about $\Phi=\pi$, changing the sign of the slope $\frac{d\mathscr{H}_{ij}}{d\Phi_{ji}}=-\frac{d\mathscr{H}_{ij}}{d\Phi_{ij}}$.. We model oscillator interactions as occurring entirely through diffusive transport of the inhibitory species bromine; oscillator $i$ is delayed by oscillator $j$ for nearly all $\Phi_{ij}$, with the the largest delay occurring just before oscillator $i$ is about to undergo a transition from the reduced to oxidized state. Experimental conditions are chosen to be consistent with the model Li et al. (2014, 2015).

State Diagram. We observed three dynamical attractors as a function of the number of arms in a star graph and the coupling strength (Fig. 3). The states are (i) unlocked (movie S5); unsynchronized oscillations of all nodes, (ii) locked (movies S2-S4); arm-hubs synchronized in-phase with a $k$-dependent phase shift between the arm and central hub, and (iii) center silent (movie S6); a non-oscillating, or intermittently oscillating central hub and unsynchronized oscillations of all arm-hubs. Using photochemical inhibition, we are also able to change state dynamically by manipulating topology. Shining strong light on a node inhibits oscillation, effectively pruning that node from the network. We were able to induce a transition from center silent to locked in a 5-arm star by shining light on two arms, effectively transforming the network to a 3-arm star while leaving the coupling constant unchanged. We also induced a transition from locked to unlocked by shining light on the hub of a 3-arm star, as shown in movies S7-S8 and Fig. S2 BZS .

In Fig. 3, experiment and theory are compared. When the coupling strength is low (large drop separation), the unlocked state is observed. A network is considered unlocked if a steady-state locking angle is not achieved during the experiment, for theory, we examine a time window commensurate with experiments corresponding to $\sim$20 oscillations. For moderate star degree, as one increases the coupling strength, phase locking is observed. Further increases to the coupling strength results in center-silent dynamics. For large star degree unlocking proceeds directly to center-silent. Conversely, as star degree is lowered, the coupling strength range over which phase-locking occurs broadens. An analytic expression for the diffusive rate $\mu~{}[\mathrm{s}^{-1}]$ is presented in the appendix, we expect $\mu$ to be inversely proportional to the drop separation and to the volume of the receiving drop BZS .

In order to match the theory-predicted state diagram with the experiment we introduced normally distributed variations of 7.5% in H+ which created a distribution in oscillator frequency. We ran each parameter combination ($k$ and $\mu$) twenty times, resampling the variation in chemistry and initial conditions each time. We also introduced a coupling strength reduction (fudge) factor $f=0.15$ to modify the predicted $\mu$. Heterogeneity is needed to both destroy synchrony when coupling is weak Strogatz (2000) and move the transition from locked to center-silent to lower $k$. Phase diagrams showing theoretical predictions without these modifications are shown in the supplement, Fig. S4 BZS .

Variations in parameters are justified because the inner, aqueous phase is composed of two reactant streams; fluctuations in flow rates can therefore produce droplets of varying composition and size. For simplicity, we consider chemical heterogeneity alone. The reduction of the effective coupling coefficient has been noted previously Tompkins et al. (2014). These results indicate our idealized model is qualitatively, but not quantitatively, correct and further, that system heterogeneity is an important parameter.

Steady State Locking Angles and Phase Model. The phase interaction function $\mathscr{H}_{ij}$, assuming only bromine transport between oscillators, is nearly completely negative meaning that coupled oscillators mutually delay one-another, Fig. 2 (bold curve). This produces a steady state phase shift between two oscillators, at the maximum phase difference of $\pi$, Fig. 2 ($k=1$ curve). Within the phase-locked regime, we experimentally observe that as the number of arms increases, the delay between when the hub and arm nodes oxidize decreases, see color-coded time traces in Fig. 1d & e and Fig. 4. We examine in detail this dependence using the numerically constructed phase model.

Since all solutions exhibit arm-locked dynamics, a consequence of the arms forming a single orbit of the star network’s graph Pecora et al. (2014), we examine the dynamics along the arm-synchronized manifold $\Phi_{1}=\Phi_{2}=\dots=\Phi_{k}$ and reduce the system to a single degree of freedom $\Phi_{0a}=\Phi_{0}-\Phi_{a}$ with arm-synchronized cluster $a$ and hub [math]. The dynamics of this quotient network are described by

[TABLE]

Examining arm-synchronized dynamics of $N$ oscillators is equivalent to considering a heterogeneous pair: a node containing a single unit of volume $V$ and a large node of volume $kV$ diffusively coupled through a conduit with $k$-times the cross sectional area $A$ of the connections in the original network, inset of Fig. 4c. This volume difference manifests as an asymmetric coupling constant $\mu_{0a}=k\mu_{a0}$. We report the system dynamics for various $k$ rescaled by the coupling strength and star degree $\mathscr{G}/\left(k\mu\right)$ so that as $k$ becomes large the plotted amplitude of the dynamics does not grow as well, Fig. 2B.

The system’s locking angles are readily given by the fixed points $\mathscr{G}\left(\Phi_{0a}^{*}\right)=0$ with stability determined by the slope $\mathscr{G}^{\prime}$: if $\mathscr{G}^{\prime}<0$, the point is stable; if $\mathscr{G}^{\prime}>0$, unstable. The results show that for two equally sized nodes, there are four fixed points, two of which are stable $\Phi_{0a}^{*}=0,\pi$. However, the basin of attraction associated with $\Phi_{0a}^{*}=0$ is so small, it isn’t visible in Fig. 2B. In contrast, the basin for $\Phi_{0a}^{*}=\pi$ is larger and deeper and therefore more accessible and robust against differences in the intrinsic frequencies of the wells, as seen experimentally Toiya et al. (2010); Delgado et al. (2011).

The phase model predicts that as the volume ratio increases from unity, the hub oscillator is more delayed by the collective action of the arm cluster, shifting $\Phi_{0a}^{*}$, Fig. 4. For example, for $k=3$, theory predicts the locking angle is reduced to $0.76\pi$. As the volume ratio increases, the fixed point continues to move until the volume ratio reaches 1:33.8 where a saddle-node bifurcation eliminates the attractor, leaving $\Phi_{0a}^{*}=0$ as the only attractor, shown as a black line in Fig. 4.

These predictions compare favorably to experiment, and the full chemical model without heterogeneity, Fig. 4. In the limit of weak coupling strength, assumed by the creation of the phase model, the locking angle is $\mu$-independent. As coupling strength is increased in the point model, predictions diverge with $k$ more rapidly. Additionally, for strong coupling ($\mu=0.1-0.5$ s*-1*) the system transitions to center-silent while for weaker coupling ($\mu=0.01-0.03$ s*-1*) the system instead transitions to in-phase synchrony as predicted by the phase model. Heterogeneity introduced into the model Eqn. 1 to reproduce the phase boundaries in Fig. 3 does not change the overall dependence of locking angle on star degree, illustrated with squares in Fig. 4.

In the limit of weak coupling, $\Phi_{0a}^{*}-\pi\propto-\log{k}$, as shown in Fig. 4. While physically $k\in\mathbb{N}$ for our system, more generally, $k$ represents the volume ratio between two diffusively-coupled oscillators and can take on any positive value. We anticipate logarithmic scaling because the deviation from $\pi$ should change sign, but not magnitude upon relabelling the quotient graph, which is equivalent to inverting the volume ratio $\log{(k)}=-\log{(1/k)}$. To identify the pre-factor, we perform a regular perturbation expansion of the locking angle with $\log\left(k\right)$ as the expansion variable, with details provided in the supplement BZS . The expansion approximates the location of the attractor of Eqn. 2 according to the scaling law

[TABLE]

which compares well to the phase model up to $\log{k}\sim\mathscr{O}\left(1\right)$. The first term, $\Phi_{0a}^{*\left(0\right)}=\pi$, arises from symmetry and is agnostic to oscillator chemistry and network topology. In contrast, the next order correction encodes information specific to both. Information about the chemical reactions enters through the oscillator’s phase response curve and coupling to adjacent oscillators through $\mathscr{H}$, while information about network topology enters through $\log{k}$. We note that while networks of repulsively coupled oscillators have been modeled simply as negatively coupled Kuramoto oscillators $\dot{\Phi}_{i}=\sum_{j=1}^{N}\mathcal{A}_{ij}\sin{\left(\Phi_{j}-\Phi_{i}\right)}$ Acebrón et al. (2005); Giver et al. (2011), our result shows explicitly that because $\mathscr{H}\left(\pi\right)=\sin\left(\pi\right)=0$, sine coupling does not predict the observed topology-dependence for branching networks.

Discussion and Conclusion. In this Letter we show how the coupling strength and topology of an inhibitor-coupled BZ reaction-diffusion star-network controls transitions between distinct dynamic states. We demonstrate that the high-dimensional reaction-diffusion system of the star network, consisting of $(k+1)\times M$ variables (Eqn. 1), can be reduced to a low-dimensional phase model of $k$ variables and further simplified to a quotient graph of only 1 phase variable that remains able to make quantitative predictions of the locking angle, Eqn. 3. The symmetry-based argument in reducing the star network to a quotient graph (Eqn. 2) requires perfectly identical oscillators. The strong agreement with experiment indicates that topology-dependent dynamics are robust against experimental imperfections. One of the ongoing objectives of the study of networks is to elucidate the connection between topology and dynamics. Our theoretical prediction of the locking angle, Eqn. 3, explicitly separates topological effects from chemical dynamics. The theory is universal and therefore readily transferable to other oscillator networks. These results demonstrate the utility of model experimental reaction-diffusion systems for both testing theories of network dynamics and providing engineering principles for dynamic soft materials.