Oregonator Scaling Motivated by the Showalter–Noyes Limit

Article pubs.acs.org/JPCA

Oregonator Scaling Motivated by the Showalter−Noyes Limit Harold M. Hastings,*,†,‡ Richard J. Field,§ Sabrina G. Sobel,∥ and David Guralnick∥ †

Division of Science, Mathematics and Computing, Bard College at Simon’s Rock, 84 Alford Road, Great Barrington, Massachusetts 01230, United States ‡ Department of Physics and Astronomy and ∥Department of Chemistry, Hofstra University, Hempstead, New York 11549, United States § Department of Chemistry, University of Montana, Missoula, Montana 59812, United States ABSTRACT: The Belousov−Zhabotinsky (BZ) reaction is the prototype oscillatory chemical system. We investigate here a new scaling of the Oregonator model of BZ chemical kinetics and use this scaling to elucidate fundamental properties of BZ dynamics. In particular, the Showalter−Noyes criterion for oscillation, that the product [BrO3−][H+] exceeds a critical value, arises naturally as a subcritical Hopf bifurcation in this setting, as does the reduction to a two-variable model. We thus provide chemical explanations of the role of time scales in the BZ reaction.

■

INTRODUCTION

dx /dt = k 3[BrO3−][H+]2 y − k 2[H+]xy + k5[BrO3−][H+]x − 2k4x 2

The prototype oscillatory chemical system, the Belousov− Zhabotinsky (BZ) reaction,1−10 is the reversible-metal-ionredox-couple-catalyzed oxidation of a readily brominated organic substrate, e.g., malonic acid, 2,4 pentanedione, or 1,4cyclohexanedione, by the bromate ion (BrO3−) in a strongly acidic (∼0.1−1 M H2SO4), aqueous medium. In a thin layer of unstirred, ferroin-catalyzed reagent one typically sees target patterns of outgoing waves of oxidation (blue) moving in a red (reduced) medium. We investigate here a new scaling of the Oregonator10 dynamical equations resulting from the Field− Körős−Noyes (FKN)2 skeleton model of BZ chemical kinetics and use this scaling to elucidate fundamental properties of BZ dynamics. The three variables in the Oregonator, denoted x, y, and z, respectively correspond to [HBrO2], [Br−], and 2[M(n+1)+], the oxidized form of the metal-ion catalyst, e.g., [Ce4+] or [Fe(phen)33+]. Three additional species, [MA], [H+], and [BrO3−], are present in significant excess, presumed to not change significantly from their respective initial values and used as parameters. Here [MA] denotes the sum of concentrations of all forms of malonic acid (including brominated and dibrominated malonic acid). The stoichiometric factor f is used to approximate the dynamics of bromide production from the oxidation of brominated and dibrominated malonic acid by M(n+1)+ and its brominated daughter products. The resulting Oregonator model10 is given by eq 1 © XXXX American Chemical Society

dy/dt = −k 3[BrO3−][H+]2 y − k 2[H+]xy + f /2kc[MA]z dz /dt = 2k5[BrO3−][H+]x − kc[MA]z (1)

These equations are typically scaled form

3,4,11,12

in the following

ε(dx /dτ ) = qy − xy + x(1 − x) ε′(dy/dτ ) = −qy − xy + f dz /dτ = x − z

(2)

The large difference between the time scales of Br− and HBrO2, given by the ratio ε′/ε ≈ 10−3, suggests adiabatically eliminating y as a variable by taking the limit ε′ → 0 in eq 2 and setting y = fz/(q + x). This yields the standard two-variable reduction of Oregonator3,4,11,12 dynamics, eq 3, which is convenient for use in large-scale calculations.13 ε(dx /dτ ) = x(1 − x) + fz((q − x)/(x + q)) dz /dt = x − z

(3)

However, the relation between the underlying BZ dynamics and the scaling and consequent model reduction of eqs 2 and 3 is not readily apparent. We therefore propose another scaling based directly upon BZ chemistry. Received: June 21, 2016 Revised: September 30, 2016

A

DOI: 10.1021/acs.jpca.6b06285 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

■

THEORY Motivation: Relationships among Sets of Chemical Species. Inspection of eq 1 shows that the parameter [BrO3−] always appears as the product p = [BrO3−][H+]. Furthermore, Showalter and Noyes14 observed that the appearance of oscillations in the BZ reaction requires that the product of concentrations p = [BrO3−][H+] must exceed a critical value, pc (∼0.045 M2 for the ferroin-catalyzed reaction in typical reagent mixtures), at least for sufficiently high [H+]. We call this critical concentration product the Showalter−Noyes limit. These observations suggest that the quantity p = [BrO3−][H+] has a fundamental significance. Similarly, the variable y  [Br−] always appears as the product u = [H+]y. Other Models. These properties of Oregonator reaction rates, that [BrO3−] always appears in the product p = [BrO3−][H+] and y  [Br−] always appears in the product u = [H+]y, hold for many other mass-action models of the BZ reaction.15−19 Although not explicit in reaction F3 involving [BrO3−] in the Showalter−Noyes−Bar-Eli model,15 this reaction is reaction 7 in the Gyorgyi−Rempe−Field model,17 where the reaction rate contains [BrO3−] in a product with [H+]. Rescaling. We shall rescale the Oregonator so that p appears only in the ratio of time scales3,4,11,12 (1/ε) of the “fast activator/autocatalytic species” (HBrO2) to the “slow inhibitor species” (M+(n+1)). This rescaling will provide a natural link between chemical and mathematical descriptions of dynamics. Rescaling x = pX, u = pU, z = p2Z, or, equivalently

ε ∼ 1/p ε′ ∼ 1/([H+]p) ε′/ε ∼ 1/[H+]

Thus, the ε′ → 0 limit of eq 2, namely, eq 3, corresponds directly with the high-[H+]/low-pH limit of eq 5 under constant p, namely, eq 7, and the scaling parameter 1/p appears naturally as the inverse of the Showalter−Noyes product. Figure 1 illustrates this limiting process.

Figure 1. High-[H+]/low-pH limit yields the reduced two-variable model: illustration. Rate constants from ref 3 (k2 = 3 × 106 M−3 s−1, k3 = 2 M−3 s−1, k4 = 3 × 103 M−1 s−1, k5 = 42 M−1 s−1), f = 1.4, parameters [MA] = 0.05M, [H+] as noted below, p = [Bro3−][H+] = 0.04 M2. Numerical solutions with the stiff equation solver in winpp.20 Initial conditions: X0 = 0.005, Z = 0.001, U at steady state with respect to X0 and Z0; thus, U0 = 3.33286 × 10−9. Graph shows oscillations in X for [H+] = 0.3 (purple), 1 (blue), 3 (green), and 10 M (red). (The rescaled HBrO2 concentration X = (1/p)[HBrO2] has units 1/M.) Note convergence to oscillations in the two-variable model (black); initial oscillations are essentially the same in all cases. This limiting process holds for both X and Z for all parameter values.

X = (1/p)x U = (1/p)u = (1/[BrO3−])y Z = (1/p2 )z

(8)

(4)

yields the following rescaled Oregonator with steady state independent of p and [H+] dX /dt = p(k 3U − k 2UX + k5X − 2k4X2)

Center Manifold Reduction. The equations in eq 3 have been interpreted as a center manifold reduction of BZ dynamics.21 This center manifold reduction to two variables appears as a natural chemical limit ([H+] → ∞, with p held constant) in the rescaling given by eqs 4−6), potentially clarifying its range of validity. Figure 2 illustrates how well the center manifold approximation Uss tracks the actual dynamics of the (adiabatically eliminated) rescaled bromide concentration U. The center manifold approximation Uss (eq 6) closely tracks the actual value of U except at the sharp rise at the onset of oscillations, where it rises more slowly (Figure 2). Figure 3 provides a further comparison of three-variable and two-variable dynamics. These models were used to further investigate dynamics near the Showalter−Noyes limit. We found that the Showalter− Noyes product p = [BrO3−][H+] serves as a bifurcation parameter for a Hopf bifurcation as p increases through a critical value critical value pc. The empirical critical values pc of p are pc ≈ 0.012 M2 in the three-variable model at [H+] = 1 M and pc ≈ 0.01 M2 in the reduced two-variable model, both less than the experimental empirical Showalter−Noyes limit of 0.045 M2 for the ferroin-catalyzed reaction.14 We always found regions of bistability for p near pc, implying a subcritical Hopf bifurcation as discussed below.

dU /dt = [H+]p( −k 3U − k 2UX + (f /2)kc[MA]Z) dZ /dt = 2k5X − kc[MA]Z

(5)

The parameters p and [H+] determine the relative time scales of the variables X, U, Z; increasing [H+] simply increases the relative speed of the “very fast” variable U. Taking the high-[H+]/low-pH limit of eq 5 with p held constant then adiabatically eliminates the variable U  (1/ [BrO3−])y, whose quasi-steady-state value is given by Uss = (f /2)kc[MA]Z /(k 3 + k 2X )

(6)

analogous to the standard reduction to two variables. This limit yields a two-variable reduction directly analogous to eq 3 dX /dt = p((f /2)kc[MA](k 3 − k 2X )Z /(k 3 + k 2X ) + k5X − 2k4X2) dZ /dt = 2k5X − kc[MA]Z (7)

More precisely, comparing eqs 2 and 5 yields the following relationships among scaling parameters B

DOI: 10.1021/acs.jpca.6b06285 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Typical results are shown in Figure 3. As shown in Figure 3A and 3B, there is a region of bistability for p near pc ≈ 0.012 M2 for the illustrated parameter values. Similar bistability occurs near the critical value of p for all parameter values; see the Discussion. Figure 3C illustrates the role of the relative time scales of (rescaled) bromide and bromous acid dynamics, ε′/ε ∼ 1/[H+], eq 8, namely, that the empirical critical value p0 is reduced as the separation of time scales approaches the high[H+]/low-pH limit. Finally, Figure 3D illustrates how the high[H+]/low-pH limit approaches two-variable dynamics, see also Figure 2 above.

■

DISCUSSION AND CONCLUSIONS

We have developed a new, simple rescaling of the Oregonator model, closely related with the actual chemistry of the underlying FKN skeleton model.2 This rescaling arises from a search for sets of chemical species which occur in all model reactions, in this case, the pair [BrO3−] and [H+] and also [Br−] and [H+]. This rescaling provides a natural link between chemical and mathematical descriptions of Oregonator dynamics. This approach may also be applied to many other mass action models for the Belousov−Zhabotinsky dynamics.

Figure 2. Illustration of the center manifold approximation. Time series of the rescaled bromide concentration generated by eq 4 U (green curve), and the corresponding quasi-steady-state projection onto the center manifold (eq 6) Uss = ( f/2)kC[MA]Z/(k3) + k2X) (black). (The rescaled Br− concentration U = [Br−]/[BrO3−] is dimensionless, as is its quasi-steady-state value Uss.) Here [H+] = 1 M; other parameters and initial values are as in Figure 1. In general, for all parameter values, the actual, very fast bromide dynamics closely tracks its quasi-steady-state value except at the sharp rise at the onset of oscillations, where it rises more slowly, reaching ∼2/3 of the quasisteady-state value in this illustrative case.

(1) Adiabatic elimination of very fast Br− dynamics (relative to HBrO2 dynamics) yielding a rescaled two-variable

Figure 3. Dynamics of the (rescaled) Oregonator (eq 4) near the Showalter−Noyes limit for oscillation. Rate constants from ref 3, f = 2, parameters [MA] = 0.05M, [H+] and p = [BrO3−][H+] as noted below, numerical solutions with the stiff equation solver in winpp,20 initial concentrations at steady state X0 = 1.99924 × 10−6 M−1, U0 = 2.09980 × 10−5, Z0 = 3.35872 × 10−3 M−3 unless otherwise noted below. (A) Bistability (associated with a subcritical Hopf bifurcation, see the Discussion) in the three-variable Oregonator near the Showalter−Noyes limit for these parameters (here [H+] = 1 M, p = 0.012 M2.14 (Black curve) Approach to steady state from X0 = 5 × 10−6 M−1. (Red and blue curves) Approach to limit cycle from X0 = 5.5 × 10−6 and 6 × 10−6 M−1, respectively. (B) Closeup view of the steady-state region in A. (C) Effect of [H+] on the Showalter−Noyes limit; here p = 0.01 M2. At [H+] = 10 M, the black curve approaches the limit cycle from X0 near steady state (to 6 significant figures) as does the red curve from X0 = 6 × 10−6 M−1. At [H+] = 1 M, the blue curve approaches steady state from same X0. (D) Comparison of the rescaled 2-variable Oregonator (heavier black curve, eq 6, the high-[H+]/low-pH limit of the full 3-variable dynamics) with the full 3-variable Oregonator (lighter black curve). Same initial conditions as black curve in C, but note that [H+] is no longer a parameter in the rescaled two-variable Oregonator. The reduced two-variable model essentially reproduces three-variable dynamics with high [H+]. C

DOI: 10.1021/acs.jpca.6b06285 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A Oregonator3,4,11,12 is equivalent to taking the limit [H+] → ∞ with p constant. (2) The concentration product p may be considered as a secondary bifurcation parameter. The Showalter−Noyes limit14 appears naturally as a critical value pc of the bifurcation parameter p and thus the ratio of activator to inhibitor time scales. The value of pc depends upon [H+], decreasing slightly and approaching a finite constant as [H+] → ∞ and also upon the stoichiometric factor f in Oregonator dynamics. (3) The decrease in the critical value pc in the two-variable model compared with the three-variable model might be related mathematically to a more rapid rate of rise of U = [Br−]/[BrO3−] at the start of each excitation in the twovariable model (Figure 1, above), making oscillations more likely given the role of U in the dynamics of the activator X (eq 5). The underlying chemistry provides an analogous explanation: the two-variable model is the high-[H+] limit of three-variable dynamics under constant p  [BrO3−][H+]. In the high-[H+] limit, reactions with mass action rates involving more moles of H+ than moles of BrO3−, in particular, reactions involving Br−, are sped up relative to other reactions, making oscillations more likely given the role of Br− in the dynamics of the activator HBrO2 (eq 1). (4) Moreover, as described below, examination of a closely related model, the generic Boissonade−De Kepper twovariable model,25,26 appears to confirm that the Hopf bifurcation as p increases through the critical value pc is subcritical, as we found in all investigations of the Showalter limit, as well as in investigations of the related FitzHugh−Nagumo neural model.27,28 More generally, eq 5 may be considered as a relaxation oscillator with fast activator variable X and slow inhibitor variable Z,22−24 which reproduces the dynamics of eq 4 at high[H+]/low-pH. In this interpretation, the parameter p controls the ratio of the activator and inhibitor time scales. Relaxation oscillations require a suitable separation of time scales, which occurs here when the concentration product p = [BrO3−][H+] exceeds a critical value pc, the Showalter−Noyes limit. This role of time scales is formalized in the generic Boissonade−De Kepper two-variable model25,26

biochemical excitable media. There are however some important differences in the stability of wavefronts in spatially extended versions of eqs 7 (BZ dynamics in unstirred media) and 9, similar to highly simplified cardiac dynamics.30 This approach to rescaling might be applied to the FitzHugh−Nagumo and related neural models.31−33 A similar explicit linkage between chemical and mathematical descriptions of dynamics can help establish ranges of validity of model reductions as well as guide experimental probes of underlying dynamics in a variety of biochemical settings.

■

Corresponding Author

*E-mail [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was partially supported by a Camille and Henry Dreyfus Senior Scientist Mentorship to H.M.H.

■

(9)

The relationship becomes clear after rescaling eq 9 as follows dx /dt = τ( −(x 3 − μ0 x + λ) − ky) dy/dt = x − y

(10)

Here x is the fast activator (our X) and y is the slow inhibitor (our Z), and the relationship between the time scale τ and Showalter−Noyes product p = [BrO3−][H+] is given by τ ∼ p((f /2)kc[MA]) ∼ p

REFERENCES

(1) Field, R. J. Chaos in the Belousov−Zhabotinsky Reaction. Mod. Phys. Lett. B 2015, 29, 1530015. (2) Field, R. J.; Körös, E.; Noyes, R. M. Oscillations in Chemical Systems. II. Thorough Analysis of Temporal Oscillation in the Bromate-Cerium-Malonic Acid System. J. Am. Chem. Soc. 1972, 94, 8649−8664. (3) Scott, S. K. Oscillations, Waves and Chaos in Chemical Kinetics; Oxford University Press: New York, 1995. (4) Epstein, I. R.; Pojman, J. A. An Introduction to Nonlinear Chemical Dynamics; Oxford University Press: New York, 1998. (5) Kiprijanov, K. S. Chaos and Beauty in a Beaker: The Early History of the Belousov-Zhabotinsky Reaction. Ann. Phys. (Berlin, Ger.) 2016, 528, 233−237. (6) Zaikin, A. N.; Zhabotinsky, A. M. Concentration Wave Propagation in Two-Dimensional Liquid-Phase Self-Oscillating System. Nature 1970, 225, 535−537. (7) Hastings, H. M.; Field, R. J.; Sobel, S. G. Microscopic Fluctuations and Pattern Formation in a Supercritical Oscillatory Chemical System. J. Chem. Phys. 2003, 119, 3291−3296. (8) Field, R. J.; Försterling, H.-D. On the Oxybromine Chemistry Rate Constants with Cerium Ions in the Field- Kö rö s-Noyes Mechanism of the Belousov-Zhabotinskii Reaction: The Equilibrium HBrO2 + BrO3− + H+ ↔ 2BrO2• + H2O. J. Phys. Chem. 1986, 90, 5400−5407. (9) Hegedus, L.; Wittmann, M.; Noszticzius, Z.; Yan, S.; Sirimungkala, A.; Försterling, H.-D.; Field, R. J. HPLC Analysis of Complete BZ Systems. Evolution of the Chemical Composition in Cerium and Ferroin Catalysed Batch Oscillators: Experiments and Model Calculations. Faraday Discuss. 2002, 120, 21−38. (10) Field, R. J.; Noyes, R. M. Oscillations in Chemical Systems. IV. Limit Cycle Behavior in a Model of a Real Chemical Reaction. J. Chem. Phys. 1974, 60, 1877−1884. (11) Tyson, J. J. Scaling and Reducing the Field-Koros-Noyes Mechanism of the Belousov-Zhabotinskii Reaction. J. Phys. Chem. 1982, 86, 3006−3012. (12) Field, R. J. Oregonator. Scholarpedia 2007, 2, 1386. (13) Merkin, J. H. On wave trains arising in the two-variable Oregonator model for the BZ reaction. IMA J. Appl. Math. 2013, 78, 513−536. (14) Showalter, K.; Noyes, R. M. Oscillations in Chemical Systems. 15. Deliberate Generation of Trigger Waves of Chemical Reactivity. J. Am. Chem. Soc. 1976, 98, 3730−3731.

dx /dt = −(x 3 − μ0 x + λ) − ky dy/dt = (x − y)/τ

AUTHOR INFORMATION

(11)

Thus, eq 9 can be considered as a polynomial approximation to the center manifold reduction (eqs 6 and 7). Note that the equations in eq 9 are in fact a special case of FitzHugh−Nagumo neural model.27,28 Glass and Mackey29 used this relationship to help study cardiac dynamics before the advent of modern supercomputers and developed a theory of D

DOI: 10.1021/acs.jpca.6b06285 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A (15) Showalter, K.; Noyes, R. M.; Bar-Eli, K. A modified Oregonator model exhibiting complicated limit cycle behavior in a flow system. J. Chem. Phys. 1978, 69, 2514−2524. (16) Györgyi, L.; Turányi, T.; Field, R. J. Mechanistic details of the oscillatory Belousov-Zhabotinskii reaction. J. Phys. Chem. 1990, 94, 7162−7170. (17) Györgyi, L.; Rempe, S. L.; Field, R. J. A novel model for the simulation of chaos in low-flow-rate CSTR experiments with the Belousov-Zhabotinskii reaction: a chemical mechanism for twofrequency oscillations. J. Phys. Chem. 1991, 95, 3159−3165. (18) Györgyi, L.; Field, R. J. Simple models of deterministic chaos in the Belousov-Zhabotinskii reaction. J. Phys. Chem. 1991, 95, 6594− 6602. (19) Györgyi, L.; Field, R. J. A three-variable model of deterministic chaos in the Belousov−Zhabotinsky reaction. Nature 1992, 355, 808− 810. (20) Ermentrout, B. README for WPP. http://www.math.pitt.edu/ ~bard/classes/wppdoc/readme.htm (accessed Aug 1, 2015). (21) Fernandez, A.; Sinanoğlu, O. Subordination of the Fast-Relaxing Degree of Freedom in the Center Manifold of the BelousovZhabotinsky System. Phys. Rev. A: At., Mol., Opt. Phys. 1985, 31, 2736−2737. (22) Sagués, F.; Epstein, I. R. Nonlinear Chemical Dynamics. Dalton Trans. 2003, 7, 1201−1217. (23) Ginoux, J.-M.; Letellier, C. Van der Pol and the History of Relaxation Oscillations: Toward the Emergence of a Concept. Chaos 2012, 22, 023120. (24) Kinoshita, S. Introduction to Nonequilibrium Phenomena. In Pattern Formations and Oscillatory Phenomena; Kinoshita, S., Ed.; Elsevier: Amsterdam, 2013; pp 1−20. (25) Boissonade, J.; De Kepper, P. Transitions from Bistability to Limit Cycle Oscillations. Theoretical Analysis and Experimental Evidence in an Open Chemical System. J. Phys. Chem. 1980, 84, 501−506. (26) Epstein, I. R.; Showalter, K. Nonlinear Chemical Dynamics: Oscillations, Patterns, and Chaos. J. Phys. Chem. 1996, 100, 13132− 13147. (27) FitzHugh, R. Impulses and Physiological States in Theoretical Models of Nerve Membrane. Biophys. J. 1961, 1, 445−466. (28) Nagumo, J.; Arimoto, S.; Yoshizawa, S. An Active Pulse Transmission Line Simulating Nerve Axon. Proc. IRE 1962, 50, 2061− 2070. (29) Glass, L.; Mackey, M. C. From Clocks to Chaos: The Rhythms of Life; Princeton University Press: Princeton, NJ, 1988. (30) Showalter, K. Quadratic and Cubic Reaction−Diffusion Fronts. Nonlinear Sci. Today 1995, 4, 1−10. (31) Izhikevitch, E. M. Dynamical Systems in Neuroscience: The Geometry of Excitability and Bursting; MIT Press: Cambridge, MA, 2007. (32) Ermentrout, G. B.; Terman, D. H. Mathematical Foundations of Neuroscience; Springer Science & Business Media: New York, 2010. (33) Gerstner, W.; Kistler, W. M.; Naud, R.; Paninski, L. Neuronal Dynamics: From Single Neurons to Networks and Models of Cognition; Cambridge University Press: Cambridge, U.K., 2014; http:// neuronaldynamics.epfl.ch/online/index.html (accessed Aug 24, 2016).

E

DOI: 10.1021/acs.jpca.6b06285 J. Phys. Chem. A XXXX, XXX, XXX−XXX