740
J . Phys. Chem. 1989, 93,140-149
shown in Figure 3, the value of kOw obtained from the formation curves was seen to increase linearly with the initial concentration of [3’,2,3] up to the maximum concentration used (6.4 X lo4 M). The absorption of the parent compound in the 400-nm region prevented the use of more concentrated solutions. In Figure 3 are also shown the values of kobd which were obtained from a computer simulation of the reaction mechanism represented by Scheme I. For these calculations, the value of the rate constant for the reduction of [3’,2,3] was that derived from the experiments shown in Figure 3 ( k = 1.5 X lolo M-I s-l); the curves describing the formation of [2’,2,3], via formation of [3‘,2,2], were always approximated to a single-exponential behavior. The comparison between the values sf kobd obtained experimentally and those obtained from the computer simulation of the
system as described above, indicates that kintramust have a value of ca. 5 X lo7 s-I or higher. This estimate is in agreement with calculations madeI5 using the Hush model and the available spectroscopic information which give a value of ca. lo9 s-l.
Acknowledgment. We thank Dr. Nadia Camaioni for solving some software problems and G. Gubellini and A. Monti for technical assistance. This work was partially supported by Minister0 della Pubblica Istruzione of Italy (Quota 40%). Registry NO. [3,2,3](PF6),, 94499-25-9; [3’,2,3](PF&, 1 1 1557-16-5; [3’,2,3’](PF6)6, 1 1 1557-14-3; [2,2,3], 94499-28-2; [2,2,2], 94499-29-3; [2’,2,2], 1 1 1557-22-3; [2’,2,2’], 1 1 1557-25-6; [2’,2,3], 1 1 1557-21-2; [2’,2,3’], 1 1 1557-24-5; HC02Na, 141-53-7; (CH,)3COH, 75-65-0; H’, 12385-13-6; COz’-, 14485-07-5.
Chemlcal Vortex Dynamics in the Belousov-Zhabotlnsky Reaction and in the Two-Variable Oregonator Model W. Jahnke,+ W. E. Skaggs,t and A. T. Winfree* Department of Ecology and Evolutionary Biology, University of Arizona, Tucson, Arizona 85721 (Received: May 25, 1988)
Like the Belousov-Zhabotinsky reaction, the two-variable Oregonator model supplemented by molecular diffusion supports rotating vortices of chemical activity in two dimensions (“rotors”) and in three dimensions (”scroll rings”). The dynamics of both kinds of organizing center are found computationally to include a regularly patterned “meander” in ways that depend on chemical recipe parameters. A search for comparable regularities in the laboratory was rewarded after detection and removal of a number of experimental artifacts.
Introduction Excitable media tend to organize themselves into periodic wave trains emanating from vortexlike ”rotors”. In three dimensions the vortex is typically a closed ring that moves under its own dynamical rules. Though predicted’ and observed2 many years ago, such vortex rings have not yet been subjected to quantitative laboratory investigation in any chemically excitable medium. Attempts to do so with the Belousov-Zhabotinsky reagent are still frustrated by technical difficulties, but it should at least be possible to preview such experiments c~mputationally~ and perhaps to sidestep some avoidable laboratory difficulties by forethought in the choice of recipe parameters. Accordingly, we undertook computational experiments4 using a qualitatively acceptable caricature of the chemical reactions involved in the BelousovZhabotinsky reaction. There are now several such models, none clearly superior in all respects to the others, and they are hard to compare quantitatively. Because some of the important parameters remain poorly measured, computation from such reaction-diffusion equations still remains to some degree an exercise in curve fitting. For present purposes we chose the earliest, simplest, and most studied model: the two-variable version of the Oregonator with Tyson’s “Lo” parameters (ref 5 ; fitted from diverse “room temperature” data). Any other kinetics could just as easily be plugged into our computational machinery: when the appropriate equations are finally agreed upon and parameters fmed at widely accepted values, these rough trial computations can be rerun to check the models by quantitative comparison with experiment. The two-variable Oregonator has solutions that represent propagation of a steep wave front by interplay of HBrO, diffusion Permanent address: Institut fuer Physikalischeund TheoretischeChemie, Universitaet Tuebingen, 7400 Tucbingen, Federal Republic of Germany. *Permanent address: Department of Psychology, University of Colorado, Boulder, CO 80309.
with an autocatalytic reaction that quickly generates HBrOz (using Br- as an intermediary species that remains always in equilibrium with local instantaneous HBr02). As in an analogous classical problem in evolutionary genetic^,^,^ the square of the propagation speed should be proportional to the product of the HBr02 diffusion coefficient, the rate coefficient of its autocatalysis, the [H+], and the [BrO,-] : this predicted relationship was confirmed experimentally in one version of the reaction, but with a rate coefficient whose value seems uncertain by a factor of 2.5 Behind the front the medium is refractory to immediate reexcitation; the recovery process is mediated by the metal ion catalyst. Our present purpose is a qualitative exploration of rotor behavior, both numerically and in the laboratory. The Oregonator has too many well-known deficiencies to make quantitative results very useful. It predicts a much weaker dependence of rotor period on [H’] than is o b ~ e r v e d . ~It does not predict phase advances observed after a pulse of Br-.8 It does not adequately distinguish the roles of malonic acid and of bromomalonic acid, which seem (1) Winfree, A. T. Sci. Am. 1974, 230(6), 82. (2) (a) Winfree, A. T. Faraday Symp. Chem. Soc. 1975,9,38. (b) Welsh, B.;Gomatam, J.; Burgess, A. Nature (London) 1983, 304, 611. (3) (a) Panfilov, A. V.; Pertsov, A. M. Dokl. Akad. Nauk. USSR 1984, 274, 1500 (in Russian). (b) Panfilov, A. V.; Rudenko, A. N . Physica D (Amsterdam) 1987, 28D, 215. (c) Winfree, A. T.; Winfree, E. M.; Seifert, H. Physica D (Amsterdam) 1985, 170, 109. (4) (a) Winfree, A. T.;Jahnke, W.; Skaggs,W., manuscript in preparation. (b) Jahnke, W.; Henze, C.; Winfree, A. T. Nature, in press. ( c ) Winfree, A. T.; Jahnke, W. J . Chem. Phys., in press. (5) Tyson, J. J. In Oscillations and Traveling Wwes in Chemical System; Field, R. J., Burger, M., Eds.; Wiley: New York, 1985; p 93. (6) (a) Fisher, R. A. Ann. Eugen. 1937, 7, 355. (b) Showalter, K.; Tyson, J. J. J . Chem. Educ. 1987, 64,142. (c) Murray, J. D. J . Theor. Biol. 1976, 56, 329. (7) (a) Smoes, M. L. J . Chem. Phys. 1979,71,4669. (b) Edelson, D. Int. J . Chem. Kinet. 1981, 13, 1175. (8) (a) Ruoff, P. J . Phys. Chem. 1984, 88, 2851. (b) Dolnik, M.; Schreiber, I.; Marek, M. Phys. Lett. A 1984, IOOA, 316. ( c ) Dolnik, M.; Schreiber, I.; Marek, M. Physica D (Amsrerdam) 1986, 2JD, 78.
0022-3654/89/2093-0740$01 .50/0 0 1989 American Chemical Society
The Journal of Physical Chemistry, Vol. 93, No. 2, 1989 741
Computations in Chemical Vortex Dynamics important in thefparameter. It takes no account of air exchange in thin unstirred layers of liquid such as are commonly used to study wave propagation, but such exchange can be seen to make a big differences9 It is not immediately adaptable to substitution of cerium by another catalyst,%JOand when used anyway to mimic wave propagation in ferroin-catalyzed reagent, it predicts more dramatic ferroin excursions than have been observed.” In fact, the range is independent of the concentration of catalyst used. RovinskyI2 corrects this and proposes an extra term and an adjustment of parameters (at 40 OC, however) and demonstrates that, as in every other excitable medium, the proposed interplay of reaction and diffusion supports stable rotors radiating spiral waves (often called “reverberators” in translations of Russian).
Computations Zero-Dimensional Local Excitation and One-Dimensional Waves. The equations of the two-variable Oregonator, for an unmixed, spatially extended chemical medium, are5,13J4 du/dt = [U - u2 - ~ v ( u- q ) / ( U + q)]/t + D1 V’U du/dt = u
- v + 0 2 V‘V
(1) with spatially uniform steady state (stable, f o r f > 1 2Il2) at - 1)’ + 4q(f+ 1))1/2] U, = U, = 0.5[(1 - -J f (cf+ (2) which for f not close to 1 and 0 C q 1 + 2ll2 and small q). The u-nullcline’s middle branch then crosses u = us, near Uthr = f q , so the threshold for excitation (“nearby”) is Uthr - us, = qcf- cf+ 1)/cf- 1)). The ranges of u and u during an excitation also depend on the parameter values. The kinetic equations above, with t taken arbitrarily small and diffusion neglected, suggest as an approximation to the case of propagating waves that in the most excited state u approaches 1 while u stays near the steady-state value and that in the most refractory state u is again near 0 and u approaches 1/(4n. We find computationally with finite t = O.Ol,f= 3.0, q = 0.002, and retaining diffusion, that in spirals u ranges from near 0 to about 0.75 (cf. 1 ) and u ranges from near 0.02 to about 0.1 1 (cf. 0.08). Changingfto 1.5, u ranges from near 0 to about 0.85 (cf. 1 ) and u from near 0.05 to about 0.20 (cf. 0.17). k l , ... k5 are the rate constants for the five pertinent chemical reactions that the Oregonator takes into account. Together with chemical recipe parameters, they are used to make all variables dimensionless:
+
(9) (a) Vidal, C. J. Stat. Phys. 1987,48,1017. (b) Showalter, K.; Noyes, R. M.; Turner, H. J . Am. Chem. Soc. 1979,101, 7463. (c) Agladze, K. I.; Panfilov, A. V.; Rudenko, A. N . Physicu D (Amsterdam) 1988,29D,3,409. (d) Winfree, A. T. Lact. Notes Biomath. 1974,2, 243. (IO) (a) Field, R. J.; Foersterling, H.-D. J . Phys. Chem. 1986,90,5400. (b) Nagy-Ungvarai, Z.;Mueller, S.C.; Plesser, T.; Hess, B. Naturwissenschaften 1988,75,87. (c) Pagola, A.; Ross, J.; Vidal, C., preprint, 1987. (d) Nagy-Ungvarai, Z.; Tyson, J. J.; Hess, B., preprint, 1988. (11) (a) Wood, P. M.; Ross, J. J . Chem. Phys. 1985, 82, 1924. (b) Mueller, S.C.; Plesser, T.; Hess, B. Science (Washington, D.C.)1985,230, 661. (c) Mueller, S.C.; Plesser, T.; Hess, B. Naturwissenschaften 1986,73, 165. (d) Mueller, S. C.; Plesser, T.; Hess, B. Physicu D (Amsterdam) 1987, 240,71. (e) Mueller, S. C.; Plesser, T.; Hess, B. Physica D (Amsterdam) 1987,240, 87. (12) Rovinsky, A. B. J . Phys. Chem. 1986,90,217. (13) Tyson, J. J. J . Chim. Phys. 1987,84, 1359. (14) Keener, J. P.; Tyson, J. J. Physica D (Amsterdam) 1986,21D,307.
u and u are local [HBr02] and [Fe(~hen),~’]scaled by appropriate chemical parameters (see Figure 2 caption). t is time scaled to units of l/k5[malonic acid bromomalonic acid], which, with k5 = 0.4 M-’s-I and [MA BMA] = 0.12 M, is 21 s; our concentrations are chosen to match the recipe in ref 15. x is distance scaled by choice of Dl = 1 to units of (D/k5[MA BMA])l/Z, which, with D = 1.5 X cm2/s and k5, [MA BMA] as above is 0.18 mm. Papers 5, 13, and 14 scale D1 = t, which additionally involves k3[H+][Br03-] in the unit of space. D1and D2 are the dimensionless diffusion coefficients of HBr02 and of ferroin; we estimate from their molecular weights (respectively 113 and 596) that they come in the ratio D2/D1 = 113/596)II3 = 0.6 so we set Dl = 1 and Dz = 0.6. t is k5[MA + BMA]/k,[H+] [BrOC], which, with the foregoing nominal values and k3 = 40 M-2 s-l, [H+] = 0.34 M, and [BrO,-] = 0.31 M, is 0.01 1 . Our computations use the value 0.01 as in ref 5, 13, and 14. With our scaling, a change of t produced by changed [H’] or [BrOC] does not affect the unit of time or space, but if changed by [MA + BMA], it affects both. With Tyson’s scaling, the same is true only for the unit of time; but the unit of space x equals ((k3[H+][Br0
