Design of localized spatiotemporal pH patterns by means of antagonistic chemical gradients

Abstract

Spatially localized moving and stationary pH patterns are generated in two-side-fed reaction-diffusion systems. The patterns are sandwiched between two quiescent zones and positioned by the antagonistic gradients of the reactants of the self-activatory process. Spatial bistability, spatiotemporal oscillations, and formation of stationary Turing patterns have been predicted by numerical simulations and observed in experiments performed by using different hydrogen ion autocatalytic chemical systems. The formation of stationary patterns due to long-range inhibition is promoted by a large molecular weight hydrogen ion binding polymer.

---

1 Introduction

The two most influential ideas about the origin of biological patterns are Turing's reaction-diffusion (RD) model¹ and Wolpert's processing of positional information.² The Turing model describes the autonomous pattern generation in a two variable activator–inhibitor model, where the effective diffusivity of the activator component is smaller than that of the inhibitor. Among the different RD patterns predicted by Turing, the stationary ones, which are formed by a spontaneous symmetry breaking process and characterized by their intrinsic wavelength, are the most striking. The biological relevance of this mechanism has been demonstrated by Kondo and coworkers on horizontal stripes of zebrafish³ and on the formation of regularly spaced transverse ridges of the palate.⁴ The Wolpert model implies the existence of a morphogen concentration gradient which is interpreted in a threshold-dependent way. This mechanism which results in distinct regions with sharp borders is widely accepted as a paradigm in developmental biology,⁵ and it can be observed in a synthetic biochemical system.⁶ Chemical RD systems are ready to form a large variety of biologically relevant spatiotemporal phenomena,^7,8 and autonomous living-like soft structures as well.⁹ The coupling of reactions which are capable of producing pH patterns with responsive materials is a quite promising approach for developments along this latter line.¹⁰

Here, we present a method to generate localized pH patterns in the presence of antagonistic gradients maintained by the counter-diffusion of different chemicals. A Wolpert type mechanism, the interplay of gradients and the self-activatory process, is used to create a spatially organized state of the system where a localized high extent of reaction zone is separated by sharp thresholds from the outer nonreactive regions. Besides this one, a flat, unstructured spatial state also can be generated. The coexistence of two stable spatial steady states, characterized by different spatial concentration distributions, at the same boundary conditions is not trivial and called spatial bistability. We have demonstrated the existence of this nonlinear phenomenon in our system, as the spatially organized state coexists with the flat spatial state over a finite range of parameters. In the high extent of reaction zone of the spatially organized state, kinetic and diffusion driven instabilities may lead to the formation of different RD patterns. The proposed approach has been tested by numerical simulations and in experiments made in two-side-fed reactors (TSFR).¹¹ In the TSFR used here the two tanks, fed continuously with different chemicals, are connected through an agarose gel (Fig. 1a). The chemicals are distributed in a way to avoid any significant reactions in the tanks. In the gel, diffusion establishes concentration gradients of the chemical species and the RD patterns develop perpendicular to that. TSFRs have been already used to study Turing patterns in the chlorite-iodide-malonic acid (CIMA) reaction,¹² and waves in the Belousov–Zhabotinsky (BZ) reaction.¹³ In both cases the input feed chemicals are the precursors of the main dynamical species, thus the patterns form in a more or less thick layer of the medium in which the simultaneous development of different spatial instabilities is possible at different positions.^14,15 This aspect makes difficult the interpretation and the control of the developing patterns. In the systems we use here, the input feed chemicals are directly involved in the positive and the negative feedback processes. This allows to create a localized zone with a single dynamic state and provides an unambiguous control of the dynamics through the variation of the boundary conditions.

Fig. 1 Sketch of the open two-side-fed gel reactor (TSFR) with a thin gel strip (a). Typical concentration profiles at the F state (b) and at the M state (c) of the gel content obtained by numerical simulations. Numerical (d and e) and experimental (f, g and h) snapshots of the stationary spatial states. The parameters used in the simulations are: a_A0 = 1, a_hA0 = 0, b_A0 = 1.5, a_B0 = 0.4, a_hB0 = 0.6, b_B0 = 0, c_X0 = 0 and w̃ = 1.0. The input feed concentrations in the experiments were: [bromocresol green]_X0 = 0.21 mM, [Na₂SO₃]_X0 = 80 mM, [BrO₃⁻]_A0 = 200 mM, [Fe(CN)₆⁴⁻]_X0 = 0 mM, [H₂SO₄]_B0 = 6 mM (f and g) and [Na₂SO₃]_X0 = 89 mM, [IO₃⁻]_A0 = 150 mM, [(NH₂)₂CS]_X0 = 5 mM, [H₂SO₄]_B0 = 7 mM (h). The pseudocolor images were derived from grayscale ones by mapping each intensity value to a color according to the presented lookup table. High intensity (yellow) color indicates a pH below 4.8.

We use Landolt-type pH oscillators, as these proved to be beneficial in previously tested RD pattern design concepts.^16–20 These systems can be described by the following reactions:²¹

A⁻ + H⁺ ⇌ HA (R1)

C + B + H⁺ → Q (R3)

S⁻ + H⁺ ⇌ HS (R4)

In this model B stands for the oxidant (e.g. BrO₃⁻, IO₃⁻ or H₂O₂), A⁻ is the main reductant (e.g. SO₃²⁻), C denotes a second reductant (e.g. [Fe(CN)₆]⁴⁻ or (NH₂)₂CS), and P and Q are products. Reactions (R1) and (R2) implement the H⁺-autocatalytic oxidation of the substrate, that is the positive feedback, while the delayed consumption of H⁺ ions (R3) and the protonation of the substrate (R1) are responsible together for the negative feedback in the chemical kinetics. Reaction (R4) represents the protonation equilibrium of polyacrylate, which is used to create long range inhibition by slowing down the effective diffusivity of hydrogen ions. Despite its simplicity this model successfully describes the general characteristics of the different spatial states and the dynamics of the systems investigated here. The aim of our work is to produce spatial bistability, spatiotemporal oscillations and Turing patterns in the presence of antagonistic gradients of certain reactants. Experimental evidence is provided by using various types of pH oscillators,²² namely the bromate-sulfite-ferrocyanide (BSF), the iodate-sulfite-ferrocyanide (FIS), the iodate-sulfite-thiourea (TuIS) and the hydrogen peroxide-sulfite-ferrocyanide (HPSF) reactions. As these reactions show sustained homogeneous oscillations only in CSTRs (continuously fed and stirred tank reactors) in the presence of nonstop inflow of the initial reactants, the appearance of RD patterns, especially Turing-type patterns, is not trivial.

2 Methods

2.1 Experimental

The experiments were performed in a TSFR, where a thin strip of a 2% agarose (Sigma-Aldrich, A2929) gel is in contact with two stirred tanks. The width (w), the length and the height of the gel strip was 2.8 mm, 23 mm and 1 mm, respectively, and the residence time in the tanks was fixed to 120 s. The experiments in the BSF, FIS and TuIS reactions were made at T = 35 °C and at T = 25 °C in the HPSF reaction. The feed solutions of the major chemicals were stored in three separated reservoirs but entered premixed into the tanks. Reactants were distributed as follows: Reservoir 1: oxidant (NaBrO₃ (Sigma-Aldrich, ≥99%) or H₂O₂ (VWR, 30% aqueous solution) or KIO₃ (Sigma-Aldrich, ≥98%)); Reservoir 2: Na₂SO₃ (Sigma-Aldrich), Na₄Fe(CN)₆·10H₂O (Sigma-Aldrich) or (NH₂)₂CS (Sigma-Aldrich) and the indicator; Reservoir 3: H₂SO₄ (diluted from 1.0 mol L⁻¹ standard solution (VWR)). Sodium-polyacrylate (Sigma-Aldrich) was added to Reservoir 2. All solutions were prepared with ion-exchanged water and chemicals were used without further purification. The input feed concentrations ([X]₀) of the reagents are indicated in the text. Here [X]₀ denotes the concentration that species X would have in the total inlet flow prior to any reaction. Tank A is fed from Reservoir 1 and 2, while Tank B is fed from Reservoir 2 and 3. At the applied pH and residence time, the content of both tanks should be considered as nonreactive, in spite of the fact that the oxidant and reductant are both fed from Tank A. In Tank B, the role of H₂SO₄ is to protonate a part of SO₃²⁻ to HSO₃⁻. To visualize the pH patterns bromocresol green (BCG) (Sigma-Aldrich, 90%, pK_a = 4.8) or bromocresol purple (BCP) (Sigma-Aldrich, indicator grade, pK_a = 6.3) was selected. The pictures were taken by an AVT Stingray F-033B (14 bit) camera, and recorded by the Streampix (Norpix) software. The camera was fixed above the experimental setup, which was enlightened by a LED source through a 530 nm long-pass filter to provide better contrast and to protect the ferrocyanide from light. We used the ImageJ program for image processing.²³

2.2 Computational details

The simulations were made by the dimensionless equations derived from the model represented by reactions (R1)–(R4). The details of the derivation are presented in the ESI.† Here, dot, ∂_t and Δ denote the dimensionless time derivatives and Laplacian. The state of the content of Tank A and B are calculated by the following dimensionless equations:

ȧ_X = −κ₁a_Xh_X + κ₋₁a_hX + a_X0 − a_X (1)

ċ_X = −κ₃b_Xc_Xh_X + c_X0 − c_X (5)

ṡ_X =−κ₄s_Xh_X + κ₋₄(s_tot − s_X) − s_X (6)

Here X refers either to Tank A or B. The variables a_X, a_hX, h_X, b_X, c_X and s_X denote the dimensionless concentrations of species A⁻, HA, H⁺, B, C and S⁻, respectively (for the definitions see ESI†). The typical dimensionless input feed concentrations in Tank A are a_A0 = 1, a_hA0 = 0 and b_A0 > 1, while in Tank B these are a_B0 = 1 − a_hB0, 0 < a_hB0 < 1 and b_B0 = 0. The input feed of C and S⁻ is the same for both tanks as c_A0 = c_B0 ≥ 0 and s_A0 = s_B0 ≥ 0. The values of κ₂, , κ₃, κ₄ and κ₋₄ are set to 5 × 10⁵, 5 × 10¹, 5 × 10³, 5 × 10¹⁰, 5 × 10⁷ respectively. The composition of the gel content is calculated by the following equations:

∂_ta = −κ₁ah + κ₋₁a_h + Δa (7)

∂_tc = −κ₃bch + Δc (11)

∂_ts = −κ₄sh + κ₋₄(s_tot − s) + 0.01Δs (12)

The dimensionless variables a, a_h, h, b, c and s correspond to the concentrations of A⁻, HA, H⁺, B, C and S⁻ in the gel, respectively. These variables are time and space dependent. We assume that the diffusion coefficient of the H⁺ in the gel is four times larger,²⁴ while the diffusion coefficient of low mobility polyacrylate is 100 times smaller than that of the other species. The concentrations at x = 0 and at x = w are the same as that of Tank A and B, respectively.

The partial differential equations were discretized with standard second order finite difference scheme on a 100×400 mesh. The resulting systems were solved by the SUNDIALS CVODE²⁵ solver using backward differentiation formula method. The absolute and the relative error tolerances were 10⁻¹² and 10⁻⁷, respectively.

3 Results and discussion

To achieve suitable conditions for the development of localized patterns we apply a uniaxial bipolar design with the antagonistic gradients of the reactants of reaction (R2), the main step of the positive feedback process. This means that B and HA are only present in Tank A and B, respectively. This distribution establishes the localization of the autocatalytic zone in the direction of the concentration gradients, because reaction (R2) proceeds only in the middle zone of the gel. The separation of these two chemicals is necessary and sufficient for our purpose, namely to generate a three-zoned structured spatial state, so the other input chemicals (A⁻, C, polycarboxylate, and the pH indicator) are fed from both sides. In the experiments the concentration of HSO₃⁻ in Tank B is controlled by the input feed concentration of H₂SO₄. This configuration results in the appearance of two different stationary states. At the one we call F (“Flow”) state, there is no significant chemical reaction inside the gel, therefore all gradients are smooth (Fig. 1b and d). In the experiments, the indicator shows its high pH color all along the gel at this state (Fig. 1f). The other spatial state is characterized by three distinct zones (Fig. 1c and e) and we designate it as the M (“Mixed”) state of the gel. The composition of the three zones is essentially different. At the left zone the concentration of A⁻ and B is high, but the concentration of HA is low. In the middle zone A⁻ is totally consumed and accordingly H⁺ appears in a relatively high concentration. The right zone is characterized by the presence of HA and A⁻ and the lack of B. The position of the maximum of the H⁺ concentration distribution (x̃*) is mainly determined by the concentration gradient of B and can be calculated as:Here, b_A0 is the dimensionless concentration of B in Tank A. This formula, which is based on the 1 : 1 stoichiometry of reaction (R2) and on a linear approximation of the gradients, is in good agreement with the numerical simulations. The width of the high H⁺ concentration band is directly proportional to w (Fig. S1 in ESI†). This type of size invariance of the three-zoned-pattern is the consequence of the bipolarity of the system.² In the experiments, at the M state of the system a sharp zone forms inside the gel, where the indicator switches to its low pH color (Fig. 1g). When iodate was used as the oxidant, iodine precipitate formation was observed within the acidic zone due to the Dushman reaction (Fig. 1h).²⁶ According to the simulations and also the experiments made in four different systems (BSF, FIS, TuIS, and HPSF), the F and the M states may coexist, that is spatial bistability. In chemical RD systems this behavior have been demonstrated only in the presence of unidirectional concentration gradients (e.g. in one-side-fed reactors), previously.^27–29

The width of the low pH zone of the M state and the parameter domain of spatial bistability decreases when the negative feedback is enhanced by the increase of the input feed concentration of C (Fig. 2a). Spatial bistability vanishes at a critical value and spatiotemporal oscillations appear with nonzero amplitude and frequency, which is signature of a subcritical Hopf bifurcation. The simulated and the experimental phase diagrams (Fig. 2a and b) show a cross-shaped topology, that is typical for open activatory–inhibitory systems.³⁰ In the oscillatory domain there are still three distinct zones in the gel, but now the central one is periodic as it can be seen in the time-space plots in Fig. 2c and d (simulation and experiment). In the presented cases the acidic signal appears and disappears abruptly. The formation of localized spatiotemporal oscillations have been observed in experiments made by bromate (BSF) and iodate (FIS and TuIS) oxidants, but not by hydrogen peroxide (see the movies in ESI†). The period of the spatiotemporal oscillations in the investigated systems is about 12–20 min at the applied conditions (Fig. 2d). The correspondence between the simulated and experimentally recorded dynamics (Fig. 2) reveals the relevance of the simple mechanism presented by reactions (R1)–(R4). Representative time-space plots of the oscillations and the local time evolution at the center of the middle zone in the FIS and the TuIS systems are shown in Fig. 3. We assume that the failure of finding oscillations in the HPSF case is connected to the rather complex kinetics of the reaction between hydrogen peroxide and ferrocyanide,³¹ which provides the negative feedback. The rate of the oxidation steps, reaction (R2) and (R3) are quite different in case of bromate, iodate and hydrogen peroxide. We have tested this effect by numerical simulations and we have found a square root scaling (Fig. S2 in ESI†) between the relevant kinetic parameters (the rate constants of reaction (R2) and (R3)) and the parameter which determines the diffusive gradients, that is the thickness of the gel. Besides, the large amplitude oscillations presented in Fig. 2 the simulations predict also small amplitude ones (Fig. S3 in ESI†), where the pH variation of the middle zone is so tiny, that the experimental observation is unlikely by using standard pH indicators.

Fig. 2 Simulated (a) and experimental (b) nonequilibrium phase diagram and time-space plots of simulated (c) and experimental (d) oscillations. The dimensionless parameters used in the simulations are: a_A0 = 1, b_A0 = 1.5 and w̃ = 2.0. The numerical time-space plot is made at c_X0 = 1.5, a_B0 = 0.4 and a_hB0 = 0.6. The experiments were made with BSF reaction, and the input feed concentrations are: [Na₂SO₃]_X0 = 80 mM, [BrO₃⁻]_A0 = 200 mM, [bromocresol green]_X0 = 0.21 mM. The time-space plot was recorded at [Fe(CN)₆⁴⁻]_X0 = 10 mM, [H₂SO₄]_B0 = 8 mM.

Fig. 3 Sustained localized oscillations in the FIS (a and b) and in the TuIS reaction (c and d). Time-space plots along the axis x (a and c) and the local oscillations at x = x* on arbitrary grayscale (b and d). Input feed concentrations: [Na₂SO₃]_X0 = 89 mM, [IO₃⁻]_A0 = 150 mM, [BCG]_X0 = 0.21 mM. In the FIS reaction [Fe(CN)₆⁴⁻]_X0 = 20 mM and [H₂SO₄]_B0 = 5.5 mM. In the TuIS reaction: [(NH₂)₂CS]_X0 = 5 mM and [H₂SO₄]_B0 = 7.0 mM.

At this point, we could adapt the design method which was originally developed to produce stationary patterns in one-side-fed reactors.¹⁶ The key points of this method are finding spatial bistability and spatiotemporal oscillations in a controlled way and ensuring long-range inhibition by using an appropriate large molecular weight complexing agent. In hydrogen ion self-activated RD systems sodium polyacrylate (NaPAA) have been efficiently used for the second purpose.³² In the model the NaPAA only reacts with hydrogen ions in reaction (R4) and its diffusion coefficient is set to be a hundred times smaller than that of the other components. According to the simulations, above a critical input feed concentration of the polycarboxylate compound, symmetry breaking patterns, a row of spots develop (Fig. 4a). The asymmetric shape of these spots is the consequence of the antagonistic gradients. In the experiments, different pattern formation scenarios could be detected above a critical concentration of NaPAA. A single row of low pH spots may appear through an oscillatory destabilization of the M state (Fig. 4b) or by step-by-step formation from the F state (Fig. 4c). In the first scenario, starting from a stable M state (Fig. 4b) the decrease of the concentration of HSO₃⁻ in Tank B leads to transient spatiotemporal oscillations. The oscillations stop as low pH spots form and are stabilized in the gel. In the second case (Fig. 4c), the spots appear in the originally F state of the gel, when the concentration of HSO₃⁻ in Tank B is increased. The wavelength of the experimentally observed pattern is about 2–3 mm.

Fig. 4 Simulated (a) and experimental (b and c) formation of Turing patterns. In each case the two snapshots show the starting and the final state, while the time-space plot represents the transitions. The dimensionless parameters used in the simulations are: a_A0 = 1, b_A0 = 1.5, c_X0 = 2.0, a_B0 = 0.4, a_hB0 = 0.6, s_X0 = 0.5 and w = 2.0. The experiments were made with the input feed concentrations: [Na₂SO₃]_X0 = 80 mM, [BrO₃⁻]_A0 = 200 mM, [bromocresol green]_X0 = 0.21 mM, [–COO⁻]_X0 = 16 mM, [Fe(CN)₆⁴⁻]_X0 = 10 mM, [H₂SO₄]_B0 = 7.4 mM (b) and [Fe(CN)₆⁴⁻]_X0 = 5 mM, [H₂SO₄]_B0 = 7.0 mM.

4 Conclusions

In this report, we presented a method to produce different types of positional chemical patterns by the interplay of appropriately designed antagonistic gradients of chemicals and the activatory–inhibitory kinetics. This arrangement could be relevant to biological systems where uniform spatial distribution is rare. Here we used it to generate pH patterns by different reactions which do not show oscillations in batch reactor only in the presence of continuous feed. We anticipate that this method can be applied to any chemical or biochemical RD system, where the diffusive gradients of the reagents of the positive feedback process can be maintained. This is the case in most of the reactions which are known as CSTR oscillators.⁷ The use of TSFRs broadens the scope of the investigation of patterns in spatially distributed systems compared to the popular one-side-fed reactors, as the separation of the reactants avoids the development of oscillations or other instabilities in the feed tanks.³³ In a pure Wolpert mechanism, the process which generates the gradients is separated from the one which interprets the positional information. In our design method, the chemicals which provide the gradients are also involved in the signaling mechanism, which makes a direct coupling between these two ingredients and results in the formation of spatial bistability. This arrangement also allows the coupling of an additional physico-chemical process to this core mechanism to induce an additional localized chemical signal. The formation of the stripe of iodine in the FIS experiments (Fig. 1f) is an example of it. We have shown that by introducing a proper chemical negative feedback process localized time-periodic pattern can be induced. The experimental confirmations were made by using four different variants of Landolt type pH oscillators. To initiate Turing instability we have added a large molecular weight hydrogen ion binding agent, in this wise we have created a system, where the Wolpert mechanism is used to localize Turing patterns. The spatially periodic positional pH pattern is an example of the rational combination of the two developmental concepts.³⁴

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

The authors thank the support of the National Research, Development and Innovation Fund (119360).

Notes and references

1. A. M. Turing, Philos. Trans. R. Soc., B, 1952, 237, 37–72.
2. L. Wolpert, J. Theor. Biol., 1969, 25, 1–47.
3. S. Kondo and T. Miura, Science, 2010, 329, 1616–1620.
4. A. D. Economou, A. Ohazama, T. Porntaveetus, P. T. Sharpe, S. Kondo, M. A. Basson, A. Gritli-Linde, M. T. Cobourne and J. B. Green, Nat. Genet., 2012, 44, 348.
5. G. Forgacs and S. A. Newman, Biological physics of the developing embryo, Cambridge University Press, 2005.
6. A. S. Zadorin, Y. Rondelez, G. Gines, V. Dilhas, G. Urtel, A. Zambrano, J.-C. Galas and A. Estevez-Torres, Nat. Chem., 2017, 9, 990.
7. I. R. Epstein and J. A. Pojman, An Introduction to Nonlinear Chemical Dynamics: Oscillations, Waves, Patterns, and Chaos, Oxford University Press New York, 1998.
8. A. S. Mikhailov and G. Ertl, Chemical Complexity: Self-Organization Processes in Molecular Systems, Springer, 2017.
9. M. Onoda, T. Ueki, R. Tamate, M. Shibayama and R. Yoshida, Nat. Commun., 2017, 8, 15862.
10. J. Horváth, Chem. Commun., 2017, 53, 4973–4976.
11. P. De Kepper, J. Boissonade and I. Szalai, in From sustained oscillations to stationary reaction-diffusion patterns, Springer, 2009, pp. 1–37.
12. V. Castets, E. Dulos, J. Boissonade and P. De Kepper, Phys. Rev. Lett., 1990, 64, 2953.
13. Z. Noszticzius, W. Horsthemke, W. McCormick, H. L. Swinney and W. Tam, Nature, 1987, 329, 619–620.
14. I. Lengyel, S. Kádár and I. R. Epstein, Phys. Rev. Lett., 1992, 69, 2729.
15. B. Rudovics, E. Dulos and P. De Kepper, Phys. Scr., 1996, 1996, 43.
16. J. Horváth, I. Szalai and P. De Kepper, Science, 2009, 324, 772–775.
17. J. Horváth, I. Szalai and P. De Kepper, Phys. D, 2010, 239, 776–784.
18. I. Szalai, J. Horváth, N. Takács and P. De Kepper, Phys. Chem. Chem. Phys., 2011, 13, 20228–20234.
19. H. Liu, J. A. Pojman, Y. Zhao, C. Pan, J. Zheng, L. Yuan, A. K. Horváth and Q. Gao, Phys. Chem. Chem. Phys., 2012, 14, 131–137.
20. I. Molnár and I. Szalai, J. Phys. Chem. A, 2015, 119, 9954–9961.
21. G. Rabai, ACH - Models Chem., 1998, 135, 381–392.
22. M. Orbán, K. Kurin-Csörgei and I. R. Epstein, Acc. Chem. Res., 2015, 48, 593–601.
23. C. A. Schneider, W. S. Rasband and K. W. Eliceiri, Nat. Methods, 2012, 9, 671.
24. G. Schuszter, T. Gehér-Herczegh, Á. Szűcs, Á. Tóth and D. Horváth, Phys. Chem. Chem. Phys., 2017, 19, 12136–12143.
25. A. C. Hindmarsh, P. N. Brown, K. E. Grant, S. L. Lee, R. Serban, D. E. Shumaker and C. S. Woodward, ACM Trans. Math. Softw., 2005, 31, 363–396.
26. S. Dushman, J. Phys. Chem., 1904, 8, 453–482.
27. P. Blanchedeau and J. Boissonade, Phys. Rev. Lett., 1998, 81, 5007–5010.
28. P. Borckmans, K. Benyaich and G. Dewel, Int. J. Quantum Chem., 2004, 98, 239–247.
29. J. Boissonade, P. De Kepper, F. Gauffre and I. Szalai, Chaos, 2006, 16, 037110.
30. J. Boissonade and P. De Kepper, J. Phys. Chem., 1980, 84, 501–506.
31. G. Rabai and I. R. Epstein, Inorg. Chem., 1989, 28, 732–736.
32. I. Szalai and P. De Kepper, J. Phys. Chem. A, 2008, 112, 783–786.
33. I. Szalai, J. Horváth and P. De Kepper, Chaos, 2015, 25, 064311.
34. J. B. Green and J. Sharpe, Development, 2015, 142, 1203–1211.

---

Footnote

† Electronic supplementary information (ESI) available: Derivation of the numerical model, additional figures & movies of spatiotemporal oscillations. See DOI: 10.1039/c8ra08028c

---