J. Phys. Chem. 1995,99, 12804-12808
12804
Analysis of a Model of Chlorite-Based Chaotic Chemical Oscillators Istviin Lengyel, Liiszld GyiSrgyi, and Irving R. Epstein* Chemistry Department, Brandeis University, Waltham, Massachusetts 02254-9110 Received: February 2, 1995; In Final Form: May 1, 1995@
A four-variable CSTR model proposed by RAbai and Orban to describe a number of chlorite-containing oscillatory reactions produces complex periodic oscillations and chaos. Detailed numerical simulations and two-parameter bifurcation diagrams with changing inflow concentrations and rate constants were calculated to reveal the dynamical properties of the model and to compare them with the results of experiments on systems containing chlorite- and sulfur-containing reductants such as thiosulfate, thiocyanate, and thiourea. It is suggested that the model requires more flexibility in describing the variable stoichiometry of the sulfurcontaining reductants and that it may need to account for the reactivity of chlorine dioxide. Experiments are suggested to establish the values of some of the parameters and to illuminate mechanistic uncertainties in the model.
TABLE 1: General Model for Chlorite Ion Based Chemical Oscillators9
Introduction The discovery of the first chlorite-containing chemical oscillators’.2 marked the beginning of a decade in which the number of known oscillating chemical reactions jumped by nearly an order of magnit~de.~ The family of “chlorite-driven” chemical oscillators contains dozens of reactions involving a wide variety of substrates4 Within this group of reactions, the richest subfamily consists of those oscillators that include an iodine species, generally iodide, iodate, or molecular iodine. However, the most complex dynamical behavior among the chlorite-based oscillators, including isothermal chaos in a flow reactor (CSTR), has been observed when the other reactants are sulfur compounds: thiosu1fate;f’ thiourea,’ and thiocyanate.8 No mechanism was proposed for these complex oscillators until 1993, when Rhbai and Orbh9 suggested a general model for the chlorite-based oscillators. In contrast to the case of oscillators involving chlorite and iodine species, where the feedback arises from the hydrolysis of iodine,I0-l2 Rhbai and Orbin suggest that in other chaotic chlorite oscillators the feedback process originates from the chlorite side. Their analysis shows that oscillations and chaos occur in the model for certain values of the rate constants, reactant concentrations, and flow rate. However, a detailed investigation of the range of dynamical behavior over an extended range of the parameters has not previously been performed. In this paper, we present a more detailed analysis of the RAbai-OrbAn model.
Proposed Model The model in Table 1 consists of six reactions and involves four variable concentrations, [C102-1, [HOCl], [C1202], and R, which represents a general reductant species whose oxidized form is RO. In a CSTR with inputs of C102- and R, the usual configuration for experiments, the differential equations corresponding to the model are
‘Abstract published in Advance ACS Abstracts, August 1, 1995.
reaction
+ R + H+ -.. HOC1 + RO (MI) CI02- + HOCl + H+ -.. Cl202 + H,O (M2) C1202 + R + H2O - 2 HOCl + RO (M3) HOC1 + R - RO + C1- + H+ (M4) C1202 + C102- - 2 Cl02 + CI- (M5) C1202 + H20 -. C101- + CI- + Cl02-
ratea VI
= 0.1b[C102-][R]
~2
= 104b[C10~-1[HOC11
~3
= (2 x 105)[C1202][R]
~I
2 H+
= ( 1 x IO4)[HOC1][R] = (2 x 104)’[C120~][C10z-] v6 = 10.5[c1202] v4 ~5
(M6)
All rate constants have units of M-I s-l, except kg, which is in s-l. These rate constants include a constant proton concentration. Corrected in ref 18. a
d[HOCl]ldt = v, - v2
+ 2v3 - v, - ko[HOC1]
d[C1202]ldt= v2 - v3 - v5 - v6 - ko[Cl20,]
(1)
where vl-vg are the rates of reactions Ml-M6 as defined in Table 1. Inflow concentrations are designated by a zero subscript, and ko is the reciprocal residence time. As noted above, the proposed model differs from previous mechanistic descriptions of chlorite-based oscillators in that the dynamical regulatory process is found on the chlorite side of the general chlorite-reductant oscillatory system. Model 1 is autocatalytic in HOCl. The autocatalytic process can be identified as the sum of reactions M2 and M3, in which each mole of HOCl consumed in reaction M2 leads to the production of two moles in reaction M3. The precursor of the autocatalyst is the reactive intermediate C1202. This species can react in three ways: (a) with the substrate R in reaction M3 to regenerate the autocatalyst, (b) with the reactant C102- in reaction M5 to give unreactive products, and (c) by hydrolysis in reaction M6 to give unreactive products. The autocatalyst HOCl is removed from the autocatalytic cycle by reaction M4 when both [HOCI] and [R] are high. Because the rates of reactions MI, M3, and M4 depend upon the properties of the substrate R and because a variety of substrates appear in the chlorite-driven oscillators, it is important to investigate how the values of the rate constants affect the properties of the model. We note that this four-variable model can be simplified to three variables by applying a steady state
0022-3654/95/2099-12804$09.0010 0 1995 American Chemical Society
J. Phys. Chem., Vol. 99, No. 34, 1995 12805
Chlorite-Based Chaotic Chemical Oscillators approximation to C1202. Although the dynamical behaviors of the three- and four-variable models are very similar, we study here only the more detailed version. We wish to investigate the effect of varying rate constants, and it is conceivable that changing the rate constants might carry us out of the range of validity ,of the steady state approximation. The three-variable model should be useful, though, in analytical studies with fixed rate constants. After summarizing our computational methods, we study the effect of varying the rate constants and investigate the parameter range in which complex dynamical behavior appears. Finally, we compare the transition to chaos in the model and in the experiments, address some stoichiometric problems, and suggest modifications and future experiments to clarify several details.
4.0
3.5 m
9x
-6 3.0 0
2.5 2.0
*."
'
I
.
,
-2
-3
-1
1
0
lOg(ko)
Computational Methods Because of the complexity of the dynamical phenomena to be investigated, we used several stiff numerical integrators and continuation packages to analyze the model. The LSODE, GEAR, and semi-implicit Runge-Kutta integrators gave identical results. We were not successful in calculating the line of Hopf bifurcations in the two-parameter bifurcation diagrams with the AUTO86 continuation package,13 but the CONT packageI4proved adequatefor the task. To determine the period of oscillation, we used numerical integration and in some cases we compared the result with that obtained from the continuation packages. Both approaches gave the same period, but numerical integration was much more effective computationally. In the bifurcation diagrams, regions of oscillation with different periods are shown with different gray levels and are denoted Pi where i is the period, Le., the number of extrema per cycle. Black areas on the bifurcation diagrams are "chaotic regions". Most of the area inside these regions shows chaotic oscillations, but in many cases there are windows of periodic oscillation as well. Choosing parameters in a chaotic region thus does not guarantee that chaotic behavior will result. In most cases, the existence of chaos was identified by following a period doubling sequence as we varied a parameter. In some cases, we also calculated the largest Lyapunov exponent. We shall see in some of the figures that odd and even period oscillations are separated by a narrow range of aperiodic oscillation, which is also an indication of deterministic chaos.
.
, -4
3.6 2.1 0
3.2 m
2.05
50
-9 x
2.8
I0 N
2.00 -2 4
0
Y
-2
-3
-1B
-2.0
22
lOs(ko)
-1
0
Figure 1. Two-parameter bifurcation diagrams in the [Rlo-ko (a) and [clO~-]o-k~(b) parameter planes. Fixed parameters: (a) (k,, k2, k3, k4, k5, k6, [c102-]0) = (0.1, 1 x lo4, 2 x lo5, 1 x 104, 2 x lo4, 10.5, 3.7 x and (b) kl, k2, k3, k4, ks, k6, [Ro]) = (0.1, 1 x 104, 2 x lo5, 1 x 104, 2 x 104, 10.5, 1.8 x The solid line is the line of saddle-node bifurcations, and dotted lines are the lines of Hopf bifurcations. indicates codimension two bifurcations.
+
-2.8 -2.9 -3.0
E
E -3.1 v
0
Analysis of the Model We attempted to select dynamically important pairs of reactions for two-parameter bifurcation diagrams where we map the ranges of multistability and simple and complex oscillations. The rate constant of reaction M1 has no significant effect on the dynamics unless it is much smaller than that of reaction M2. At significantly higher values of kl no autocatalysis or oscillations are observed. In the relevant parameter range, the dynamically important pairs of parameters appear to be (a) the inflow concentrations of the reactants, [Cl02-]0 and [R]o; (b) the rate constants k2 and k4 of the reactions in which, respectively, HOC1 is recycled and removed from the autocatalytic cycle; (c) the rate constants k3 and k4 of the reactions in which, respectively, R feeds and suppresses the autocatalytic cycle; (d) the rate constants k3 and k5 of the reactions in which, respectively, C1202 produces the autocatalyst and is removed from the autocatalytic cycle. Because all the variables oscillate in this model, the flow has a very important role, supplying chlorite and R in every oscillatory cycle. The reaction is a true CSTR oscillator. A natural bifurcation parameter is the reciprocal residence time, k ~ .In Figure la and b we show two-parameter bifurcation
-0
-3.2 -3.3 -3.4 -2.80 -2.78 -2.76 -2.74 -2.72 -2.70 -2.68 -2.66 -2.64 -2.62 -2.60
Wl6) Figure 2. Bifurcations of periodic orbit with changing flow rate. Fixed parameters: (kl, k2, k3, k4, k5, k6, [Cl02-]0, [Ro]) = (0.1, 1 x 104, 2 x 105, 1 x 104,2 x 104, i o . 5 , 5 x 10-3,5.05 x 10-3).
diagrams in the ko-[R]o and ko-[C102-]0 planes. Oscillation and bistability occur over wide ranges of flow rates. Within the oscillatory range (dotted line) the oscillations undergo period doubling and period adding bifurcations that lead to complex periodic and chaotic regions. In Figure 2 we show a representative bifurcation diagram for the amplitude of oscillations as ko changes. The transition to chaos occurs through period doubling bifurcations on the left side and period adding bifurcations on the right side. To obtain a two-parameter bifurcation diagram that also contains the range of oscillations with different periods, we calculated several dozen figures like Figure 2 with one of
Lengyel et al.
12806 J. Phys. Chem., Vol. 99, No. 34, 1995 1.6
__ - - -. SN ............. HB
I
1.5
FqsSYms
P1
"R!n P2
1.4
m
P3 P"
- F
(
1.3 n I0N
0 2
1.2
\
n
U
Y
1.1
1 .o
. .
.- . .
0.9
0.8 -3.0
.. ... .
!
.
!
.
O
-2.6
.
I
-2.2
.
I
.
I
.
I
.
-1.8
I
-1.4
.
I
.
I
-1.0
,
I
.
I
-0.6
log([clo;lo) Figure 3. Two-parameter bifurcation diagram in the [R]o/[ClOz-]olog([ClO~-]o)parameter plane. The parameter vector is (kl,kz, k3, k4, ks, kh, ko) = (0.1. 1 x 104, 2 x IO', 1 x lo", 2 x 104, 10.5.4 x 10-j).
the parameters fixed. The resolution of the second parameter was chosen according to the rate of change of complex behavior. From the bifurcation diagrams in Figure 1, we see that the optimal value of ko is about 4 x s-l, which is close to the reciprocal residence time used in the experiments, so ka was chosen with this or a close value for other bifurcation diagrams. In Figure 3 we show how the inflow concentrationsof chlorite and reductant affect the dynamics. For clarity we plot the ratio of these two parameters vs log [Cl02-]0, because the concentrations vary over a very wide absolute range but must be close to each other in order for complex dynamical behavior to occur. The rate constants of reactions M2 and M4 are closely related to the changes in the inflow concentrations. If we focus on the fate of the autocatalytic species HOCI, we may think of an increase in [Cloz-]~as roughly equivalent to an increase in k2. A similar correlation holds between [R]o and k4. In Figure 4 we show the effects of varying k2 and k4. Note that, like the inflow concentrations, k2 and k4 must be nearly equal to each other even though their absolute magnitudes can be changed by several orders of magnitude. Chaotic behavior is easiest to obtain if these rate constants are between 10" and lo5 M-I S-I and their ratio is very close to 1. Changing the other parameters does not significantly alter the range of these two rate constants in which complex behavior is found. The bifurcation diagram in this plane is open to the right: even with diffusion-limited rate constants (10" to lo1*M-' s-I), the bifurcation pattern is the same as at k2 = k4 = lo7 M-I s-I. The very narrow ranges of [ R ] ~ [ C ~ O ~and - ] Qk4/k2 over which complex behavior occurs suggest discrepancies between the model and experimental chlorite oscillators. With thiosulfate as the substrate, the [C102-]d[R]o ratio where chaos appears has been found to be 1.2 by Orbdn and Epstein5 and 1.6-2.3 by Maselko and Epstein.6 In the case of thiourea, this ratio is -2.0-2.5,7.1s while for thiocyanate it is -2.* Clearly, actual chaos can occur over a sizable range of reactant ratios. In the
-
model, however, this range is only about I%, independent of the values of the other parameters. Moreover, the rate constants of reactions M2 and M4 must be almost identical (Figure 5). Since reaction M2 involves chlorine species only, while reaction M4 has R as a reactant, it seems highly unlikely that kdk2 will be so close to unity for several different substrates. Peintler et a1.I6 investigated the reaction between chlorous acid and hypochlorous acid between pH 5 and 6. They fitted parameters k2 and kdk6 to almost a thousand experimental points and obtained values of (1.12 f 0.01) x lo6 M-2 S-I and (5.4 f 0.2) x 104 M-I, respectively. Rhbai and Orbdn used a rate constant of 10" M-I s-I for k2, assuming a buffered reaction medium at about pH 3.5. The chlorite-containing chaotic oscillators operate at various pHs. The chlorite-thiocyanate and chlorite-thiosulfate oscillators exhibit chaos at about pH 4, while the chlorite-thiourea reaction is chaotic at pH 2.5; consequently, the effect of changing this rate constant on the dynamical behavior is important. A low value of k2, corresponding to high pH, prevents oscillation unless k4 is decreased by a similar magnitude, but then only simple rather than complex oscillations are found (Figure 4). To increase [R]o and [Cl02-]0 to compensate for the lower rate constants is not a satisfactory solution, because their concentrations are already close to the experimental values. In Figure 5 we show the kdk3 vs log(k5) bifurcation diagram. Complex oscillations appear if 3.9 < log(k5) 5.6. At high values of kdk3 the limit cycle collides with the saddle-nodepoint and becomes unstable, leaving only a single stable steady state or a pair of bistable steady states. As ks is increased within the region of complex oscillation, we observe P1 -P2-P3-P4P5- ... period adding sequences. These sequences are always accompanied by a very narrow range of chaos. In mechanistic terms, the reason that there is so little freedom to vary the rate constants while maintaining the complex dynamical behavior is that the precursor of the autocatalyst, C1202, must be preserved in the autocatalytic cycle. The optimum value of ks in this model would be between 10" and IO5 M-I s-l, with k3 an order of magnitude larger than k ~These . choices ensure that most of the C1202 is recycled into the autocatalytic pathway rather than producing unreactive products (C102 and Cl-) in reaction M5. In Figure 6 we show the effects of k3 and k4, the rate constants of the two processes in which R feeds and kills the autocatalytic cycle. If k4 > 1.5 x 10" M-' s-l, then no oscillations occur, because the limit cycle collides with the saddle-node point and becomes unstable, leaving only steady states as attractors. The ratio of these two rate constants can be varied more freely than those of k2 and k4 or of the input concentrations [R]o and [Cl02-]0; however, the value of k3 cannot be changed significantly without eliminating the complex oscillations. We saw in Figure 5 that k4 can be increased arbitrarily at fixed k3 so long as k2 is simultaneously increased to keep kdk2 constant. Here, however, regardless of how k2 and /Q are chosen, k3 cannot be varied by more than about a factor of 3. As we increase k4 at fixed k3, we observe period doubling bifurcations leading to chaos and then complex period adding sequences. The higher the value of k4, the higher the period of the oscillations that appear immediately to the right of the line of chaotic behavior in Figure 6.
Conclusions The Rdbai-Orbdn model does represent the first successful attempt to simulate complex oscillations and chaos in sulfurcontaining chlorite-based oscillators. However, when the
J. Phys. Chem., Vol. 99, No. 34, 1995 12807
Chlorite-Based Chaotic Chemical Oscillators 1.20
1 .1 5
I
---
.............
1 .oo
-"
m
. I 0.95 0.90
2
3
5
4
7
6
lo&*) Figure 4. Two-parameter bifurcation diagram in the khlkz-log(k2) parameter plane. Fixed parameters: kl, k3, ks, k6, ko, [ClOs-]o, [Rlo) = (0.192 x 105.2 x 104, 10.5.4x 10-3,5 x 10-3.5 x 10-3). 5
---
0.16 -
1
SN
............. HB 0.15 -
7 0.14 -
0.13 0.1 2
4 -
I
~
- CI.
b
r
x
P1 P2 P3 p4 P5 P6 P7 P8 p9
cc
+
3 -
22
--
- SN
............ HB R"m P1
lm" P2 P3 P4 P5
i
0.08 .
0.07
-. . . . .
f'/-
2 -
.
CH
I-
1
3.5 3.7 3.9 4.1 4.3 4.5 4.7 4.9 5.1 5.3 5.5 5.7 5.9
'
0.8
'
1
0.9
'
1
1.0
'
1
1.1
'
'
1.2
'
1
1.3
'
1
'
1.4
1
'
1.5
1.6
Figure 5. Two-parameter bifurcation diagram in the kdk3--log(ks) parameter plane. Fixed parameters: (kl,kz, k4, k6, ko, [CIOZ-IO,[Rlo) = (0.1.I x 104, I x IOJ, 10.5.4 x 10-3.5 x 10-3, 5 x 10-3).
krXIOJ Figure 6. Two-parameter bifurcation diagram in the k r k 4 parameter plane. Fixed parameters: (kl,k2, ks, k6, ko, [C102-]0,[Rlo) = (0.1, I x 104,2 x 104, 10.5,4 x 10-3, s x 10-3. s x 10-3).
calculations described here are compared with experiment, one must conclude that further work remains to be done in both modeling and experiments. The only experimental twoparameter bifurcation diagram for these systems was determined in the [ClO2-]0-b parameter plane of the chlorite-thiosulfate reaction? In Figure lb, we see that more complex behavior occurs at high flow rate and high [Cl02-]0. As we change the flow rate, we observe period-doubling sequences in the model. In the experiments, however, complex oscillations arise with decreasing flow rate, and instead of period doubling there is period adding with narrow regions of chaos separating the
different complex periodic oscillations. The model and the experiments differ qualitatively. The model produces complex oscillations and chaos within narrow ranges of the reactant ratio, while the experiments allow much larger variation. The stoichiometry of this simple model is not sufficiently flexible if we would like to describe the oxidation of a specific substrate. It may be possible to correct this flaw by including additional steps characteristic of the individual substrates. The different reductants have significantly different reactivities, and their oxidations have variable stoichiometries that produce different reactive intermediates and
lOg(ks)
..
12808 J. Phys. Chem., Vol. 99, No. 34, 1995
products as the oxidantheductant ratio changes. For example, in the chlorite-thiosulfate reaction, tetrathionate and sulfate ions are generated in a variable ratio that depends upon the ratio of reactant concentrations and on the pH. Tetrathionate also reacts further with chlorite, hypochlorous acid, and chlorine dioxide. Dithio compounds are comparatively stable reactive intermediates in their oxidations and can provide intemal autocatalysis independent of the complexity of the oxidant side. Reactions M1, M3, and M4 constitute a necessarily oversimplified treatment of these and other stoichiometric aspects of the chlorite-based oscillators. On the chlorite side, chlorine dioxide is not simply a product of these reactions but also an intermediate. It can react with R and with other reductant side intermediates in a variety of ways. In the chlorite-thiourea reaction in a closed system, several extrema are observed in the absorbance of ClO2,’ demonstrating the intermediate character of this species. In the literature there is a discrepancy in the rate of formation of the key intermediate C1202. Peintler et a1.I6 measured a rate constant of 1.2 x lo6 at pH 5-6 while Rfibai and Beck” used a value about 2 orders of magnitude smaller to describe chlorine dioxide formation in the chlorite-iodine reaction at pH = 2. These differences suggest that there are may be a change in the mechanism of this reaction. Our preliminary experiments on chloritehypochlorous acid in acidic medium give results similar to Rfibai and Beck’s and indicate significant changes in the mechanism on going from pH 5 to 2. A detailed investigation will be published elsewhere. The idea introduced by Rfibai and collaborator^^^'^ that the key feedback process in certain chlorite-containing oscillators occurs on the chlorite rather than on the reductant side is highly original and quite promising. The fact that complex oscillations and chaos can be generated by this simple model is extremely suggestive. There are only two other systems, the Belousov-ZhabotinskyI9 and the peroxidase-oxidaseZ0 reactions, in which chemically realistic models have been proposed that describe dynamical behavior of this complexity. It is therefore worth probing further the extent to which this model can be amplified to describe the details of individual chlorite-substrate oscillators. Experiments to measure two-parameter bifurcation dia-
Lengyel et al. grams in the substrate-flow rate and chlorite-substrate subspaces should lead to a more detailed picture of the system and can provide information about rate constants and mechanistic details. Experiments to clarify the mechanism of the chloritehypochlorous acid reaction at low pH’s may help to clear up discrepancies between experiments done under different conditions and should provide a stronger grounding for this model or new variants.
Acknowledgment. This work was supported by the National Science Foundation (Grant CHE-9023294) and by US.Hungarian cooperative grants from the NSF and the Hungarian Academy of Sciences. We thank Kenneth Kustin for helpful discussions. References and Notes (1) De Kepper, P.; Epstein, I. R.; Kustin, K. J. Am. Chem. SOC. 1981, 103, 2133.
(2) Orbin, M.; De Kepper, P.; Epstein, I. R.; Kustin, K. Nature 1981, 292, 816. (3) Epstein, I. R.; Kustin, K.; De Kepper, P.; OrbBn, M. Sci. Am. 1983, 248, 112. (4) Epstein, I. R.; Orbin, M. In Oscillations and Travelling Waves in Chemical Systems; Field, R. J., Burger, M., Eds.; Wiley: New York, 1985; p 257. (5) Orbin, M.; Epstein, I. R. J. Phys. Chem. 1982, 86, 3907. (6) Maselko, J.; Epstein, I. R. J. Chem. Phys. 1983, 80, 3175. (7) Alamgir, M.; Epstein, I. R. I n f . J . Chem. Kinet. 1985, 17, 429. (8) Doona, J. C.; Doumbouya, S. I. J . Phys. Chem. 1994, 98, 513. (9) Ribai, G.; Orbin, M. J. Phys. Chem. 1993, 97, 5935. (10) Kern, D. M.; Kim, C.-H. J. Am. Chem. SOC. 1965, 87, 5309. (1 1) Epstein, I. R.; Kustin, K. J. Phys. Chem. 1985, 89, 2275. (12) Citri, 0.; Epstein, I. R. J . Phys. Chem. 1987, 91, 6034. (13) Doedel, E. J.; Kemevez, J. P. Software for Continuation and Bifurcation Problems in Ordinary Differential Equations. Technical Report, Applied Mathematics; California Institute of Technology: 1986. (14) Marek, M.; Schreiber, I. Chaotic Behavior in Deferminisfic Dissipative Systems; Cambridge University Press: Cambridge, 1991. (15) Doona, C. J.; Blittersdorf, R.; Schneider, F. W. J . Phys. Chem. 1993, 97, 7258. (16) Peintler, G.; Nagypil, I.; Epstein, I. R. J . Phys. Chem. 1990, 94, 2954. (17) Ribai, G.; Beck, M. T. Inorg. Chem. 1987, 26, 1195. (18) Ribai, G. J. Phys. Chem. 1994, 98, 5920. (19) Gyorgyi, L.; Field, R. J. Nature 1992, 355, 808. (20) Degn, H.; Olsen, L. F.; Perram, J. Ann. N.Y. Acad. Sci. 1979,316, 623. JP9503371
