Phase-frequency model of strongly pulse-coupled Belousov-Zhabotinsky oscillators

Among his many contributions to nonlinear science, Showalter has done groundbreaking work on elucidating the behavior of coupled chemical oscillators. Three decades ago, he pioneered the use of beads of ferroin-loaded cation-exchange resin as a medium for generating and studying novel spatiotemporal patterns in the Belousov-Zhabotinsky (BZ) reaction.¹ He later utilized this technique to create a chemical analog of the phenomenon known in biology as quorum sensing.² More recently, he has explored phase clusters³ and chimera states⁴ in these systems, which consist of large numbers of coupled BZ oscillators. Here, we present an approach to analyzing the behavior of coupled oscillators, drawing on the BZ reaction as our inspiration and source of experimental data.

Pulse-coupled oscillators represent a special class of coupled systems, in which the interactions between oscillators are limited to brief time periods (pulses). A well-known example is the nervous system, where neurons interact by means of short electrical pulses (action potentials) and pulsatile release of neurotransmitters.⁵ Theoretical investigations of pulse-coupled oscillators commonly employ simple phase models, in which oscillators are represented by their phases, and the interactions between oscillators are characterized by phase response curves (PRCs),⁶ which measure the phase change as a function of the phase of perturbation. This simple but powerful approach has been used to explain common synchronization phenomena such as in-phase^7,8 and out-of-phase⁹ oscillations, bursting,¹⁰ oscillator death,¹¹ and dynamical clustering in networks of neurons¹² as well as some system-specific behaviors such as alternating firing order.¹³ 

Phase models^14–17 were first employed to investigate synchronization in large sets of weakly coupled oscillators. Their applicability to pulse-coupled oscillators rests on the assumption that the coupling is sufficiently weak that each oscillator returns to its original unperturbed limit cycle after completing the cycle during which the pulse is received; i.e., the oscillator quickly “forgets” the perturbation.

When the coupling is strong, as in chemical or neural oscillators, this assumption may be violated:^18,19 the oscillatory dynamics is altered such that the unperturbed cycles are different before and after the perturbed cycle. Thus, the oscillators are not memory-less,⁶ which may bring about more complex collective behaviors.

Inhibitory coupling has been shown to produce in-phase (IP), out-of-phase (OP), and oscillatory-suppressed (OS) modes with two coupled phase oscillators.⁹ In addition to this, we observed complex temporal patterns (1:N, N:M, where P:Q signifies that each cycle contains P peaks of one oscillator and Q of the other) in an experimental system comprised of two pulse-coupled Belousov-Zhabotinsky (BZ) oscillators with asymmetric inhibitory coupling, that is, when the coupling strengths are not equal. We attributed the emergence of the complex patterns to long-lasting changes in the chemical composition (memory).²⁰ 

Here, we present an augmented phase-frequency model that can reproduce the observed long-term effect of inhibitory pulse-perturbations in a single BZ oscillator as well as the complex dynamics seen when two BZ oscillators are coupled via unequal inhibitory pulses.

A single pulse-perturbed BZ oscillator²⁰ is constructed as follows: Two reagent solutions (0.4 mol/l $NaBrO3+0.5$ mol/l $H2SO4$ and 0.1 mol/l malonic acid + 2 mmol/l ferroin + 0.2 g/l Triton X-100 surfactant) are pumped into a beaker of volume $Vr=15$ ml at room temperature (21 $°$C). The reaction mixture is stirred at 800 rpm, and the volume is held constant by removing the excess reaction mixture using a suction pump; the reactor residence time is 800 s. A Pt electrode is immersed directly into the reaction mixture, and a Ag$∣$AgCl$∣$KCl electrode is connected to it through a junction holding saturated $K2SO4$ solution. The potential (⁠$E$⁠) between electrodes is measured using a pH meter connected to a PC through a USB data acquisition board. The redox potential is monitored and recorded at a sampling rate of 2 Hz using a computer program, which also carries out the pulse-perturbations $via$ solenoid valves that release small volumes (⁠$Vp<$ 0.1 ml) of the inhibitor stock solution (⁠$[KBr]stock=10--200$ mmol/l containing 0.25 mol/l $H2SO4$⁠). The instantaneous concentration of the added KBr is calculated as $[KBr]=[KBr]stockVpVr−1$⁠.

The unperturbed ferroin-catalyzed BZ reaction exhibits strong relaxation-type oscillations with a relaxed period $T0$ of $100±5$ s and a waveform that resembles the action potentials of tonically active neurons. Upon perturbation with KBr, the next peak is delayed (inhibition) and the unperturbed cycles that follow are extended as well. The relaxed frequency (⁠$Ω=T0−1$⁠) is recovered after the oscillator remains unperturbed for several cycles [Fig. 1(a)]. We tested the effects of multiple perturbations by perturbing consecutive cycles after a fixed delay (⁠$τp$⁠) following each peak. Figure 1(b) shows that the perturbations have a cumulative effect: the higher the number of perturbations, the greater the deviation from the relaxed frequency.

We capture this behavior using a simple two variable model

We define the phase as the time elapsed since the last peak normalized to the relaxed cycle length, $ϕ=(t−t∗)/T0mod1$⁠, where $t$ and $t∗$ are the time since the start of the experiment and time of the last peak, respectively. The phase advances at the instantaneous frequency $ω$ and is reset to zero when $ϕ$ reaches 1 [$δ()$ is the Dirac delta function]; $ω$ slowly relaxes to $Ω$ with a time constant $τω≫Ω−1$⁠.

During any unperturbed time interval between peaks, $ϕ(t)$ and $ω(t)$ can be calculated by integrating Eqs. (1) and (2)

From here on, $ϕ$ and $ω$ refer to the variables and $ϕ(t)$ and $ω(t)$ to the integrated forms.

In order to obtain $τω$ from the experimental data such as that shown in Fig. 1(b), we need to derive the explicit form of the cycle length of an oscillator with $ω<Ω$⁠. We first obtain the inverse of Eq. (3); i.e., we solve for the time, $t$⁠, as a function of the phase, $ϕ$⁠,

where $W()$ is the Lambert W-function, defined by $z=W(zez)$⁠. To make our notation cleaner, we define $ϕ˚(t)$⁠, $ω˚(t)$⁠, and $t˚(ϕ)$ as the respective functions with initial parameters $ϕ0=0$ and $t0=0$⁠. We define the following function to calculate the unperturbed cycle length:

We find the period of the first cycle as $Θ(ω0)$⁠. To do the same for the second cycle, we first need to calculate the instantaneous frequency at the beginning of that cycle, $ω[Θ(ω0)]$⁠. Then, we can use Eq. (6) again. For any cycle after the first one in a relaxation sequence, the period can be calculated as a composition of functions

Since no perturbation occurs, all cycle lengths depend only on $τω$ and $ω0$ for the first cycle.

We determined the model parameter $τω$ by fitting the calculated periods to the relaxation segments in $T$ $vs.$ time series similar to those shown in Fig. 1(b). In each such experiment, $ω0$ of the first unperturbed cycle needed to be fitted as well, because $ω$ is an unobserved variable. Although we measured $Ω$ at the beginning of the experiments, the goodness of fit (measured by the F-test) was better when we fitted $Ω$ as well. We used the L-BFGS-B algorithm²¹ implemented in the R software package²² to obtain the fits.

In Table I, we report the parameters of fits that had $p<0.002$⁠. Our experimental results do not show an obvious dependence of $τω$ on $[$KBr$]$ or on the delay between peaks and perturbations (⁠$τp$⁠) preceding relaxation. The average (⁠$765±57s$⁠) of the results shown in Table I is close to the residence time of the reactor (800 s). Changing the threshold of $p$ to $10−4$ only slightly changes the value obtained. We use $τω=800s$ in our calculations.

Next, we extend our model by adding terms that describe the effect of perturbations

Perturbations occurring at $t=tp$ reset the phase according to the PRC, $g(ε,ϕ)$⁠, parameterized by the perturbation strength, $ε$⁠, which is equal to $[$KBr$]$ in our experiments. They also reduce the instantaneous frequency by $κ(ε)ω$⁠, where $κ(ε)$ is the perturbation strength-dependent frequency reduction constant.

We assume that $κ(ε)$ is independent of the phase of perturbation. Therefore, we can obtain $κ(ε)$ by fitting the relaxation segments in the PRC measurement period $vs.$ time series [Fig. 2(a)]. Each such sequence starts with $ω0=Ω$⁠. After the first perturbation $ω=(1−κ)Ω$⁠, followed by frequency relaxation until the next perturbation when it decreases further by $κω$⁠. In order to find $κ$⁠, we fitted the periods calculated using Eq. (6) to the experimental periods using the least squares method. In Fig. 2(b), we report the $κ$ values and a linear fit of $κ(ε)$ in the range of $1×10−5$M to $5×10−4$M. Note that at the highest perturbation strength, each perturbation was measured to decrease the instantaneous frequency by about 17%.

Besides the frequency change, oscillators experience phase reset when perturbed, which is calculated as $Δϕ=1−T/T0$⁠, where $T$ is the length of the perturbed cycle. The magnitude of $Δϕ$ depends on both the phase of the perturbation and its strength ([KBr]). The phase reset as a function of the observed phase is referred to as the latent PRC. We measured the latent PRCs of the ferroin-catalyzed BZ oscillator at $[KBr]=5×10−6$M to $5×10−4$M [Fig. 3(a)]. All PRC values are negative (phase delay), and the values decrease monotonically as $ϕ$ and $ε$ increase. PRCs at higher $[$KBr$]$ exhibit values that no longer obey $Δϕ≥−ϕ$⁠, which means that some perturbations reset the phase to negative values. Similar behavior was found in the Wang-Buzsáki neuron,²³ a model²⁴ commonly used in computational neuroscience.

We obtained the parameters required to calculate the relaxation dynamics; thus, it is possible to extract the true PRCs from the experimental data. We first calculate the instantaneous frequencies at the beginning of the perturbed cycles and immediately after the perturbations. The same formula can be used to calculate how much the phase advances from the beginning of the cycle until a perturbation (⁠$Δtp=tp−t∗$⁠), and from the perturbation until the end of the perturbed cycle (⁠$Δtl=tend−tp$⁠). The phase reset is then given by Eq. (11)

PRCs corrected for the deviation from the relaxed cycle are shown in Fig. 3(a). They express the phase reset when an oscillator with relaxed frequency is perturbed with KBr. These possess two major differences compared to the latent PRCs: (1) $Δϕ$ is smaller and (2) their magnitudes reach a plateau at about $ε=1.8×10−4$ as illustrated in Fig. 3(b) in the $Δϕ$$vs.$$ε$ plot at $ϕ=0.7388$⁠.

To generate a family of PRC functions over a continuous range of $ε$⁠, we used a functional form derived from the chemical kinetics of the BZ reaction²⁵ 

The constants $k$ and $ε0$ depend on the experimental control parameters (chemical composition, flow rate, etc.); $T0$ is the relaxed cycle length. In order to match the corrected experimental PRCs, we set $T0=100s$⁠, and we obtained the other two parameters by fitting Eq. (12) to the corrected PRCs in the concentration range of $1×10−5$M to $1.68×10−4$M (just before the corrected PRCs plateau). PRCs calculated using the fitted parameters are shown in Fig. 3(a). The maximum magnitude of the inhibitory PRCs is achieved at $ε=1.68×10−4$⁠; corrected PRCs at higher $ε$ are identical.

The variable frequency has an important impact on the peak delays caused by perturbations: the lower the oscillator’s frequency, the more the next peak is delayed by identical phase resets. The consumption of added KBr proceeds according to a pseudo-first order chemical reaction,²⁵ the rate constant of which depends only on the reagent concentrations of the BZ reaction, which are kept constant by the flow. Therefore, consuming the added KBr should take the same time, provided that the additions occur at the same phase. We need to ensure that this behavior is preserved in our model by rescaling the phase reset when $ω<Ω$⁠.

Using Eq. (5), we can find the time it takes for the initially relaxed oscillator to advance the amount of phase given by the PRC when the frequency is reduced to $ω0,r=[1−κ(ε)]Ω$ by the perturbation. We then calculate the phase reset that causes an equivalent peak delay when the oscillator has been perturbed earlier $ω0,p=[1−κ(ε)]ω$ using Eq. (3). The rescaled phase reset with the composition of the PRC, Eqs. (3) and (5), and the frequencies of the relaxed and the perturbed oscillator immediately after perturbation are given by

To model the behavior of two coupled oscillators, we used two pairs of the following generic equations:

where $εj→i$ is the coupling strength for the connection $j→i$ and $τ$ is an optional time delay between each peak of oscillator $j$ and the ensuing perturbation to oscillator $i$⁠.

The model equations are uncoupled between peaks and perturbations, which allows for the use of the integrated formulas during these segments. Thus, solutions can be obtained by an event-based method.

Starting from random initial conditions (⁠$ϕ1$⁠, $ϕ2$⁠) and oscillators at their natural frequencies, $ω1,2(t=0)=Ω1,2$⁠, we first use Eq. (5) to find the times (⁠$t1∗$⁠, $t2∗$⁠) it takes for each oscillator to reach $ϕ=1$⁠. The first firing occurs (⁠$t∗$⁠) at $min(t1∗,t2∗)$⁠, which sets the window of this iteration. We set the phase of the oscillator that fires first to zero and calculate the phase of the other oscillator using Eq. (3). We must also calculate the frequency of each oscillator at (⁠$t∗$⁠). Then, we can calculate the effect of the perturbation on the oscillator that did not fire. We set its phase and reduce its frequency according to Eq. (13) and $ωr=1−κ(ε)ω$⁠, respectively. This completes the first iteration, the results of which are used as initial conditions for the second iteration. This procedure allows us to find the peaks and perturbations when the coupling is instantaneous (⁠$τ=0$⁠).

When there is a delay between peaks and perturbations (⁠$τ>0$⁠), we need to maintain a list of upcoming perturbations, where each element contains the time of the perturbation, the target oscillator, and the strength of the perturbation. Thus, the previous procedure must be modified slightly: at the beginning, when the times of the next firings are calculated, we must check whether a peak or a perturbation occurs first. If a peak occurs first, we proceed according to the procedure above. If a perturbation occurs first, we set the time window to the time of the next perturbation (⁠$tp$⁠). Then, we calculate the phases and frequencies of the oscillators using Eqs. (3) and (4) at $tp$⁠. We perform the delayed perturbation on the target oscillator: we reset the phase and reduce the frequency as described above. Finally, we remove the perturbation from the perturbation array and start the next iteration with the results of the current iteration used as initial conditions.

We implemented an algorithm based on this approach in MATLAB²⁶ and calculated the time series for a large matrix of coupling strength pairs. In Fig. 4, we show domains of the simulated behaviors in the $ε1→2$ $vs.$ $ε2→1$ phase plane. Our model is able to simulate the behaviors displayed by two pulse-coupled BZ oscillators with unequal inhibitory coupling.²⁰ The calculated patterns appear in the correct coupling strength ranges when we use parameters extracted from experiments with a single pulse-perturbed oscillator. The 1:N and OS patterns predominate [Fig. 4(a)], when there is no time delay between peaks and perturbations (⁠$τ=0$⁠). When $0<τ<T0/2$⁠, the N:M type patterns are also reproduced [Fig. 4(b)].

Our approach is similar to the phase-amplitude description of neural oscillator models by Wedgwood et al.;²⁷ however, the amplitude is not a reliable observable in our system, which would hinder our attempts to extract model parameters from our experimental data by fitting. Our approach requires that the long-term dynamics of the frequency follow a simple exponential decay with a decay-constant (⁠$τω$⁠) that is independent of the history of the system. More complex long-term dynamics (e.g., multiple processes affecting the frequency, or decay of higher order) will require extensive modifications, potentially limiting the derivation and the use of closed form solutions. While our method assumes that the frequency of the oscillator is reduced after a strong perturbation, it is straightforward to modify it to capture frequency increase as a result of a perturbation.¹⁸ 

The description of dynamics using our model does not assume in-depth knowledge of the underlying chemistry of the system. The functions derived specifically from the chemistry of the BZ oscillator can be replaced by functional approximations. This method may straightforwardly be adapted to model the dynamics of systems (e.g., networks of neurons), where the coupling is strong and interactions cause the dynamics of the oscillators to undergo long-term changes. Since these equations can be solved on the basis of events, larger scale systems, which have been of interest in physics and computational neuroscience, can be easily explored.

This work was supported by the National Science Foundation (NSF) under Grant No. CHE-1362477 and by the Hungarian Academy of Sciences (Grant No. OTKA K-119360).