Reaction-Diffusion Dynamics of pH Oscillators in Oscillatory Forced Open Spatial Reactors.

Studying the effect of coupling and forcing of oscillators is a significant area of interest within nonlinear dynamics and has provided evidence of many interesting phenomena, such as synchronization, beating, oscillatory death, and phase resetting. Many studies have also reported along this line in reaction-diffusion systems, which are preferably explored experimentally by using open reactors. These reactors consist of one or two homogeneous (well-stirred) tanks, which provide the boundary conditions for a spatially distributed part. The spatiotemporal dynamics of this configuration in the presence of temporal oscillations in the homogeneous part has not been systematically investigated. This paper aims to explore numerically the effect of time-periodic boundary conditions on the dynamics of open reactors provided by autonomous and forced oscillations in the well-stirred part. A simple model of pH oscillators can produce various phenomena under these conditions, for example, superposition and modulation of spatiotemporal oscillations and forced bursting. The autonomous oscillatory boundary conditions can be generated by the same kinetic instabilities that result in spatiotemporal oscillations in the spatially distributed part. The forced oscillations are induced by sinusoidal modulation on the inflow concentration of the activator in the tank. The simulations confirmed that this type of forcing is more effective when the modulation period is longer than the residence time of the well-stirred part. The use of time-periodic boundary conditions may open a new perspective in the control and design of spatiotemporal phenomena in open one-side-fed and two-side-fed reactors.

Introduction

The concepts of forcing and coupling are central to the theory of nonlinear oscillators, and they have many practical applications. The appearance of these dynamical phenomena in nonlinear chemical systems has been studied widely since the 1980s.[1,2] Forcing, which is the limiting case of asymmetric coupling, can result in many interesting phenomena, for example, phase resetting, frequency locking, quasiperiodicity, and chaos control.[1,3−5] Coupled oscillators can be synchronized,[6] but they can produce many other behaviors such as beating and oscillator death, cluster and chimera states, and quorum sensing. These collective dynamics have been experimentally explored in an array of coupled chaotic electrochemical oscillators,[7] in diffusively coupled, nanoliter volume, aqueous drops containing the reactants of the oscillatory Belousov–Zhabotinsky (BZ) reaction,[8] and in coupled discrete chemical oscillators.[9−11]

One of the most remarkable features of nonlinear chemical systems is their capability to produce spatiotemporal patterns.[12] Since the observation of chemical waves in the BZ reaction,[13] a wide variety of spatiotemporal phenomena have been observed experimentally and described theoretically.[14] One way to create new types of patterns is to investigate novel, often biologically relevant experimental configurations. In these, numerical simulations often promote understanding the effects of the applied initial and boundary conditions.[15,16] However, periodic forcing and coupling of reaction–diffusion systems are still at the forefront of research. Frequency-locking, standing-wave patterns, and tongue-shaped regions of resonance have been reported in a light-sensitive form of the BZ reaction–diffusion system.[17,18] Localized clusters may form by using global photochemical feedback on the same reaction.[19] Periodic illumination has been successfully applied to control the formation of stationary Turing patterns[20] in the chlorine dioxide–iodine–malonic acid (CDIMA) reaction–diffusion system.[21] Interestingly, a recently published theoretical paper suggests the use of time-dependent boundary temperature to control the development of Turing patterns.[22] In coupled spatially extended systems, the development of superposition and superlattice patterns were reported in the CDIMA reaction.[23]

Remarkably, the experimental observations of forcing and coupling in chemical reaction–diffusion systems have been often made in open one-side-fed reactors (OSFRs, also called continuously fed unstirred reactors)[19,21,23] or in open two-side-fed reactors (TSFRs).[17,18] The central part of these reactors is a piece of porous material, often a hydrogel, which acts as a convection-less medium for the development of reaction–diffusion phenomena (Figure ).

Figure 1

Sketch of an open OSFR (a) and an open TSFR (b). Here, c, c(, c0, cA, and cB denote the concentrations in the input flow, in the gel, in the CSTR, and in tanks A and B, respectively.

Hereafter, we will call it as the gel part of the reactor. In an OSFR configuration (Figure a), the porous medium is in contact with the content of a continuous stirred-tank reactor (CSTR). The CSTR content feeds the porous material with fresh reagents, which allows maintaining sustained, far-from-equilibrium conditions. In a TSFR configuration, the gel is sandwiched between two separated open tanks of reagents (Figure b). Both configurations can be effectively used to explore chemical pattern formation.[24] In these reactors, the CSTR or the separated tanks provide the boundary conditions for the reaction–diffusion systems. In general, it is pretty apparent to apply time-independent, fixed boundary conditions. However, in the context of forcing of spatially distributed systems, it would be interesting to explore the effect of time-periodic boundary conditions too. There is a notable lack of studies investigating this direction either experimentally or theoretically.

The primary aim of this paper is to explore numerically the effect of time-periodic boundary conditions in an OSFR and in a TSFR. We selected a simple but chemically relevant model of pH oscillators,[25] one of the most widely studied nonlinear chemical systems,[26] especially in the context of pattern formation.[27,28] The pH oscillators are relatively simple redox reactions, which produce large-amplitude pH oscillations in a CSTR, and hydrogen or hydroxide ions play a crucial role in the positive feedback of their mechanism. Their robust pattern forming capacity has been experimentally demonstrated in OSFR and TSFR configurations and recently in a gel reactor with flow through channels.[27,29,30] We explore an autonomous and a nonautonomous way to generate time-periodic boundary conditions in an OSFR. The autonomous way is due to the capability of the reaction to produce oscillations in a CSTR. The oscillatory state of the CSTR content can be used to set time-periodic boundary conditions for the gel. The other method, the nonautonomous method, is to apply a sinusoidal perturbation on the input feed concentration of hydrogen ions, when the unperturbed CSTR content is on a stationary, low-extent-of-reaction state. In a TSFR, we apply the similar perturbation in one of the tanks.

Results and Discussion

At first, we recall the basic features of OSFR dynamics of the system at time-independent boundary conditions, which have been discussed in detail previously.[28,31] To discuss the dynamics of an OSFR, we must assign the actual state of the CSTR part and that of the gel content. To avoid any additional sources of instabilities, here, we set the diffusion coefficients of all species equal. When the input feed concentration of the activator (h) is below 0.031, the CSTR content is on a state, at which the overall extent of the reaction is low (therefore, it is referred to as the “flow” state). The reaction–diffusion system in the gel has two stationary states, which differs in their spatial concentration profiles (Figure a,b, see Supporting Information for the construction details of these graphs). The F state of the gel content is characterized with a flat profile (Figure b), while the important feature of the M (“mixed”) state is the outbreak of hydrogen ion concentration in the depth of the gel. The stability domains of these states overlap, that is, spatial bistability (Figure a). The dynamics of the reaction–diffusion system in an OSFR is strongly determined by the thickness of the gel (l). The domain of spatial bistability decreases, and spatiotemporal oscillations appear with the increase of l. It is important to notice that oscillations in the gel form at the stationary state of the CSTR content (Figure c). The typical frequency of these autonomous spatiotemporal oscillations in the gel is 2.5 at the conditions used here, meaning that the period is shorter than the residence time of the CSTR, that is, 1 in the nondimensional model. The domain of spatiotemporal oscillations overlaps with the domain of the M state. Therefore, the periodic behavior in the gel can be reached from the F state by increasing h.

Figure 2

Dynamics of the OSFR without forcing: a nonequilibrium phase diagram (a), the spatial profiles of the stationary states at h = 0.028 in the gel when the CSTR is on the flow state (b), and space-time plot of the spatiotemporal oscillations in the gel when the CSTR is on the flow state and the local time evolution of h in the CSTR (h0) and at the impermeable wall (h0.60) at h = 0.026 (c). F, M, and O denote the stationary states and the oscillatory state of the gel content, respectively. The parameters used in the simulations are c = 0.12, b = 1.5.

A natural way to create time-periodic boundary conditions is to set the value of h > 0.031, where the CSTR content starts to oscillate, at the conditions used here. The frequency of these CSTR oscillations is about 1/3 (Figure ). During these relaxation oscillations, a slow increase is followed by a sharp jump in the h0. The actual dynamics of the gel content at three representative values of l is presented in Figure . In a thin gel, l = 0.1, the CSTR oscillations are followed by the oscillations of gel content without any significant time delay (Figure a). At l = 0.2 (Figure b), the peaks of h0.2 in the depth of the gel appear slightly before the peaks of h0 in the CSTR. The width of the peaks is also more expansive in the gel than that of the CSTR content. When the value of l reaches the level at which oscillations may develop in the gel even at the stationary CSTR, the dynamics becomes more complex (Figure c). During the slow increase of h0, e.g., between t = 2 and 4 in Figure c, the gel content paths through a complex sequence starting from an F state-like composition, followed by spatiotemporal oscillations in the depth of the gel, and then reaches an M state-like composition. This sequence is ended at the sharp peak of h0 caused by the CSTR oscillations. The period of the resulting complex oscillations in the gel is determined by that of the CSTR oscillations.

Figure 3

Dynamics of the OSFR when the CSTR content oscillates: spatiotemporal oscillations driven by the CSTR oscillations at l = 0.1 (a), at l = 0.2 (b), and at l = 0.6 (c). The parameters used in the simulations are c = 0.12, b = 1.5, and h = 0.035.

Time-periodic boundary conditions can be made by applying a sinusoidal perturbation on the input feed concentrations. In the absence of any reactions, a perturbed CSTR is described by the following equation:where τ, c0 and c, A, and f are the residence time, the concentrations in the CSTR and in the input flow, the amplitude, and the frequency of the perturbation, respectively. The solution of this equation iswhere ω = 2πf and C is a constant. The amplitude of the forced concentration oscillations in the CSTR is smaller than A and decreases with the increase of ω, which means that the CSTR damps the forcing.

We have selected the concentration of the activator to perform the forcing as the constant value of h is replaced with the time-periodic one: h = h* + A sin(2πft). We kept h* below 0.031, that is, the limit of the development of autonomous oscillations in the CSTR at the applied conditions. At these conditions, the CSTR content is on a periodically perturbed F state. The amplitude of the forced oscillations in the CSTR decreases with the forcing frequency, in agreement with the above general result (Supporting Information, Figure S1).

The periodically forced CSTR sets time-periodic boundary conditions for the gel, where the reaction–diffusion system operates. In Figure , the observed complex spatiotemporal oscillations are presented at three representative values of the forcing amplitude. At low forcing, amplitude and frequency modulation of the spatiotemporal oscillations can be observed (Figure a). At the maximum of the CSTR oscillations, the frequency and the spatial amplitude of the oscillations in the gel increase, but the local concentration amplitude decreases. The opposite changes happen at the minimum of the CSTR oscillations. The increase of the forcing amplitude, presented in Figure b, results in forced bursting phenomena, where quiescent periods and active phases alternate. This dynamic behavior is well known in neuroscience.[32] A burst of spikes characterize an active phase. Forced bursting appears when a spiking (oscillatory) subsystem is slowly driven above and below the spiking threshold. The forcing amplitude and frequency can control the length of the quiescent periods and the number of spikes. In the domain of bistability between spatiotemporal oscillations and the M state of the gel, the increase of the amplitude of the forcing may induce a transition from the oscillations to a periodic M state, as is presented in Figure c. The periodic forcing stabilizes the M state at the expense of the autonomous spatiotemporal oscillations. The critical forcing amplitude at which this transition occurs increases with the forcing frequency (Supporting Information, Figure S2), in agreement with the previously mentioned damping effect of the CSTR on the forcing. The sinusoidal perturbation results in the decrease of the domain of bistability (compare Figure a and Figure S3 of Supporting Information). When the forcing frequency is close to the inverse of the residence time of the CSTR, the autonomous oscillatory state in the gel can be stable even at A = h*, which means that h oscillates between 0 and 2h*.

Figure 4

Dynamics of the OSFR when the CSTR is under sinusoidal perturbation (h = 0.026 + A sin(2πft)): modulation at A = 0.002, f = 0.2 (a), forced bursting at A = 0.003, f = 0.2 (b), and transition from the oscillatory state to the modulated M state at A = 0.004, f = 0.2 (c). The parameters used in the simulations are c = 0.12, b = 1.5, and l = 0.6.

In Figure , the observed spatiotemporal oscillations are presented at three representative values of the forcing frequency. At a forcing frequency, for example, f = 0.6 (Figure a), which is significantly lower than that of the autonomous oscillations, forced bursting can be observed. The increase of the forcing frequency, presented in Figure b, results in complex oscillations. In the presented case, a block of two spikes is followed by a block of three spikes, and this pattern repeats. The oscillations in the gel are not affected by the forcing at a frequency that exceeds that of the autonomous oscillations, for example, f = 3.0 (Figure c). We did not find modulation or synchronization between the CSTR and the gel dynamics in this case. The frequency of the oscillations in the gel is 2.5, like in the autonomous case (Figure c), while that of the forcing is 3.0. The appearance of these different types of oscillatory behaviors in the frequency–amplitude plane is shown in Supporting Information (Figure S2).

Figure 5

Dynamics of the OSFR when the CSTR is under sinusoidal perturbation (h = 0.026 + A sin(2πft)): forced bursting at A = 0.008, f = 0.6 (a), complex oscillations at A = 0.008, f = 1.0 (b), and simple oscillations at A = 0.008, f = 3.0 (c). The parameters used in the simulations are c = 0.12, b = 1.5, and l = 0.6.

Similar sinusoidal forcing can be applied in a TSFR, for example, by perturbing the input feed concentration in one of the tanks. The autonomous TSFR dynamics of pH oscillators has been explored previously both experimentally and numerically.[29] The most important features are summarized in Figure . Appropriately separated mixtures of reactants are continuously feed the two tanks. Tank A is fed by A–, H+, and C, and the boundary conditions are set at x = 0. Tank B is fed by A–, B, and C, and the boundary conditions are set at x = l. The nonequilibrium phase diagram of the gel content shows the standard cross-shaped topology:[33] the stability domain of the two spatial states F and M (Figure b) overlaps at small values of l, and above a critical value of l, spatiotemporal oscillations develop (Figure c). The spatiotemporal oscillations form in the middle region of the gel. The frequency of the autonomous oscillations in the gel is 2.5 at the applied conditions.

Figure 6

Dynamics of the TSFR without forcing: a nonequilibrium phase diagram (a), the spatial profiles of the stationary states at h = 0.85 (b), and the space-time plot of the spatiotemporal oscillations and the local time evolution of h in tank A (hA) and in the middle of the gel (h0.5) at h = 0.80 (c). F, M, and O denote the stationary states and the oscillatory state of the gel content, respectively. The parameters used in the simulations are c = 0.40, b = 1.5, l = 1.0.

We have applied a sinusoidal perturbation on the input flow of tank A as h = h* + A sin(2πft). In tank A, only reactions and R3 take place. Since H+ is bounded by A–, the extent of reaction is very small. Figure presents the increase of forcing frequency at a constant amplitude.

Figure 7

Dynamics of the TSFR when tank A is under sinusoidal perturbation (h = 0.80 + A sin(2πft)): modulation at A = 0.1, f = 0.3 (a) and A = 0.1, f = 0.5 (b) and simple oscillations at A = 0.1, f = 3.0 (c). The parameters used in the simulations are c = 0.40, b = 1.5, and l = 1.0.

The damping of the tank on the forcing amplitude is also important in this configuration. Low-frequency forcing (f = 0.3), presented in Figure a, may result in a complete periodic excursion through the different states of the gel content: starting from the F state, the increase of hA induces a transition to spatiotemporal oscillations and the M state. At a higher forcing frequency (Figure b), forced bursting phenomena can be induced, similarly as in the OSFR case (Figure a). In the TSFR configuration, we found synchronization phenomena, as presented in Figure c. When a high-frequency forcing (f = 3.0) is applied, the frequency and the phase of the spatiotemporal oscillations nicely adjust to that of the forcing. The sinusoidal perturbation of the TSFR results in the decrease of the domain of bistability (compare Figure a and Figure S4 of Supporting Information).

Conclusions

The current study aimed to investigate the impact of time-periodic boundary conditions in open spatial reactors. In the created experimentally realistic situations, we have varied and forced the input feed concentrations of the chemicals. In the case of an OSFR, time-periodic boundary conditions can be realized by autonomous oscillations in the CSTR part or by the periodic perturbation of the stationary state of the CSTR. However, in a TSFR, only the periodic forcing can be applied as the primary reagents are fed into separated tanks. Our investigation shows that this is a more flexible way to study forcing, but the damping effect of the tank on the forcing amplitude is a significant limitation. Similar results have been recently published in homogeneous systems (pH oscillators) with periodic forcing of the inflow rates of the reagents.[5] It was demonstrated that if the modulation period of the forcing was longer than the residence time of the chemical species in the CSTR, the system emulated the time period/frequency of the forcing. However, when the frequency of the forcing approached the natural frequency (in the absence of forcing) of the oscillations, the chemical oscillatory system exhibited other phenomena such as resonance and beats. When the time period of the forcing was too low compared to the natural time period of the system and the residence time of the chemical species, the forcing did not affect the characteristics of the system. In addition to this, similar behavior was found in a neutralization reaction when an acid solution was titrated periodically in antiphase with an alkaline solution in a CSTR by using various periodic inflow rate functions of the reagents.[34] Our simulations, made in spatially extended systems, confirmed that periodic boundary conditions created by low-frequency forcing on the CSTR part might result in the modulation of spatiotemporal oscillations and forced spatiotemporal bursting in the gel part. Synchronization can also be observed at high-frequency forcing in a TSFR. This study indicates that time-periodic boundary conditions in open spatial reactors may result in new reaction–diffusion dynamics. A limitation of our study is that we performed only 1D simulations to clarify the essential dynamics. Considerably more work will need to be done to determine the 2D effects of time-periodic boundary conditions, and experimental verification would also be necessary.

Model and Numerical Method

The formal chemical equations of the Rábai model of pH oscillators are the following:[25]Here, B stands for an oxidant (bromate, iodate, or hydrogen peroxide), A– is the unprotonated form of a weak acid (most often the sulfite ion), which is oxidized to an unprotonated form of a strong acid (sulfate ion), C denotes a second substrate, and P and Q are products. The corresponding rate equations are written asHere, κ1, κ–1, κ2, κ2′, and κ3 are the corresponding rate coefficients. Reaction represents the (+) feedback, whereas the (−) feedback is provided by reaction . The model is ready to simulate the temporal dynamics in a CSTR and the spatiotemporal dynamics in an OSFR or in a TSFR. We used nondimensional concentrations, where a, ah, h, b, and c correspond to A–, HA, H+, B, and C, respectively.

In an OSFR, the equations for the content of the CSTR part are the following:where a0, a, h0, b0, c0, h, b, and c are nondimensional concentrations in the CSTR and in the input feed, respectively. The reaction–diffusion equations for the content of the gel in 1D are written aswhere a, ah, h, b, and c are nondimensional concentrations in the gel. We applied Dirichlet boundary conditions or time-periodic boundary conditions at the gel/CSTR surface, for example, a( = a0, and no flux boundary conditions at the gel/impermeable wall surfaces, for example, (∂a)() = 0. The values of the nondimensional rate coefficients are κ1 = 5 × 1010, κ–1 = 5 × 105, κ2 = 5 × 105, κ2′ = 50, and κ3 = 1 × 103. We assume that the diffusion coefficients of A–, HA, B, and C are the same. The derivation of the nondimensional equations and parameters are detailed in Supporting Information. The H+ concentrations in the gel at a particular x are denoted in the text as h.

The TSFR dynamics can be described by the following dimensionless equations:

Tank A

Tank Bwhere aA, ah,A, hA, cA, aB, ah,B, hB, bB, cB, h, b, and c are nondimensional concentrations in tanks A and B and in the input feed of the tanks, respectively. The reaction–diffusion system in the gel part of a TSFR is described by eqs 8–12 with Dirichlet boundary conditions or time-periodic boundary conditions at surface A, for example, a(x = 0) = aA, and Dirichlet boundary conditions at surface B, for example, a(x = l) = aB. The H+ concentrations in the gel at a particular x are denoted in the text as h.

The partial differential equations were discretized with a standard second-order finite difference scheme on an equidistant grid having 200 gridpoints. The resulting set of the ordinary differential equations was solved by the SUNDIALS CVODE[35] solver using the backward differentiation formula method. The absolute and relative error tolerances were 10–10 and 10–5, respectively.