J. Chem. Soc., Faraday Trans. I, 1988, 84(2), 575-580 The Phase Response of the Explodator Matild Eszterle Department of Chemical Engineering, Institute of Physics, Technical University of Budapest, H-1521 Budapest, Hungary Zoltan Noszticziust Department of Physics, The University of Texas at Austin, Austin, Texas 78712, U.S.A. Zoltan A. Schelly" Department of Chemistry, The University of Texas at Arlington, Arlington, Texas 7601 9-0065, U. S. A. A re-examination of the 'limited Explodator' shows that with the appropriate numerical values of the parameters, this model can reproduce the experimentally observed phase response behaviour of the Belou- sov-Zhabotinskii reaction perturbed by the addition of Ag+ ion. It is shown that the phase-response method is not a crucial diagnostic test for the acceptability of a model.The response of systems to perturbation is routinely used to analyse their dynamics in many fields of sciences and engineering. In chemical kinetics, relaxation methods and flash photolysis represent well known examples. The general methodology is applicable to both linear and non-linear systems. In non-linear chemical systems such as oscillating reactions exhibiting limit cycle behaviour, the idea of examining system response is attractive since any clue obtainable in addition to usual kinetic information is greatly needed for the elucidation of the mechanism of these complex reactions. In recent papers1V2 the phase-response technique was used to test two skeleton mechanisms proposed for the Belousov-Zhabotinskii3 (BZ) oscillating reaction.The essence of the technique is the perturbation of the limit cycle system by the sudden change of the concentration of a key species and the measurement of the ensuing phase shift as compared with the unperturbed system. Periodicities of concentrations of limit- cycle systems are rather stable against perturbation in amplitude and frequency, but not in p h a ~ e . ~ Thus, in principle, comparison of the experimentally measured phase shift with that computed for a given model may be used to test the suitability of the model. In the BZ reaction, the experimental phase shift can conveniently be determined by monitoring the appearance of the next concentration maximum of the Ce4+ ion following the perturbation.2 Depending on whether the Ce4+ peak appears sooner or later than in the unperturbed system, one registers a negative or positive phase shift, respectively.The period of the unperturbed limit cycle is usually used as the unit for the phase shift and the phase of the perturbation. Ruoff compared' his experimental phase response diagrams2 (phase shift us. phase of perturbation) with those simulated for the Oregonator5 and the Explodator' models of the BZ reaction. In the experiments with which we are concerned, the perturbation of the limit cycle was achieved by the sudden addition of Ag+ ion, causing an instantaneous drop of the Br- concentration which led to the phase shift. A characteristic feature of the resulting phase response diagram is a linear portion of the curve (the 'excitable t On leave of absence from the Technical University of Budapest.575576 Phase Response of the Explodator branch ')2 with unit slope and negative phase shift. Based on the ability of the Oregonator to reproduce qualitatively the excitable branch, and the inability of the Explodator to do so with a given6 set of parameters, Ruoff concludedl that the Explodator is unable to describe essential features such as excitability'* and Ag+-induced oscillations9 of the BZ reaction. In this paper we re-examine the phase response of the Explodator and shall show that : (i) with the appropriate limiting reactions and parameters the Explodator is able to reproduce the experimentally observed excitable branch and (ii) without the knowledge of the actual rate constants it is not justified to rely on the shape of the phase-response diagram in deciding about the quality of a model, since the simulated phase response is obviously also a function of the numerical values of the parameters of the model used.Re-examination of the Limited Explodator Structural Features The Limited Explodator model6 consists of a generalized Lotka-Volterra type core (the (E 1 ) Explodator) X+Y+Z (E 2) Z + (1 +BY (B = 0.5) (E 3) y + {B} (E 4) {A}+X + (1 + a ) X (a = 1) and one or more limiting reactions, for example or 2X-,iZ+{A} {A) + X. It has two consecutive autocatalytic reactions in a serial network structure,lO in contrast to the one autocatalytic step and parallel network of the Oregonator. Noszticzius et al. (NFS) suggested6 a possible chemical identification of the Explodator variables as X 3 HBrO,, Y: Br,, Z E 3HOBr when applied to the BZ reaction.Phase Response of the Limited Explodator The dynamics of the Explodator core with the two limiting reactions (L 1 ) and ( L 2) is described by the following system of differential equations using dimensionless parameters (D) I i = F+x-xy-kx2 p = 1.5az-by-xy i = - az + xy + kx2/6. The case k = 0 corresponds to the situation where only (L 2) is considered as a limiting reaction and (L 1) is disregarded, as was done by NFS' to illustrate the limit-cycle behaviour of the model by a numerical example (not a simulation). For this purpose the following numerical parameters were chosen: a = 100, b = 1, F = 0.01, k = 0. The same parameters were also used in previous simulations of the phase response of the Explodator mode1.l The perturbation can be represented by the following 'instantaneous ' reaction Ag+ + Br, + H,O + AgBr + HOBr + H+ (1)M .Eszterle, Z . Noszticzius and Z . A . Schelly 577 that in terms of the variables of the Explodator is Ag++ Y -+ iZ. (EP) Now let us introduce the following symbols: [Ag+], for the virtual concentration of Ag+ after the perturbation, but prior to the beginning of reaction (EP) [this quantity is the same as the ‘added unreacted Ag+ ions’ in ref. (2)]; ryl-, for the Br, concentration at the instant of the perturbation but prior to the reaction with Ag+ in (EP); w]+, for the Br, concentration immediately after the ‘ instantaneous ’ reaction (EP). Since the duration of the perturbation and reaction (EP) is considered zero, clearly, the function [Y] has a discontinuity at the time t = 0 of the perturbation.[XI-,, [XI+,, [ZIP,, [Z],, and [Ag+]+, are quantities analogous to those defined above. Since X is unaffected by the perturbation, [XI-, = [XI+,. For the other concentrations, the relationship between their values prior to and after the perturbation is determined by the stoichiometry of reaction (1) and the amount of Ag+ used in the perturbation. Three qualitatively different cases are possible : (1) [YIP, > [Ag+],, i.e. small perturbation where the concentration of Br, exceeds that of the added Ag+. In this case the relationships P I + , = P I - 0 - [Ag+l, [ZI+o = P I - 0 + $[Ag+I, [Ag+l+, = 0 hold and one simply has to solve the system of differential equations (D) with the initial values of [XI+,, [v+, and [Z],, and the given parameters (a = 100, b = 1, F = 0.01, k = 0). Since Ce4+ does not appear in the model, the phase shift for component X was computed in ref.(1) and also here. (2) w]-, < [Ag+], < 2M-, + 3[Z]-,, i.e. large perturbation where the added Ag+ instantaneously reduces the Br, concentration essentially to zero (as determined by the solubiIity product of AgBr). If the size of the perturbation is in this interval, the Br, produced from Z in reaction (E3) can ultimately use up the excess Ag+. Nevertheless, so long as fy]+O E 0, the original system of differential equations with the given6.’ parameters (D’) I A? = 0.01 + x - x y j = 1502-y-xy i = -1ooz+xy must be replaced by R = 0.01 + x j = O 2 = - lOOz + (1 50/3)2 = - ~ O Z (since the 1.5 Y formed from 1 Z are instantaneously converted to Z by the Ag+).After the perturbation (2) and from this initial value [Z] decreases according to i = -502. Also, after the perturbation (3) From this intial value, the excess Ag+ concentration is used up by Y as fast as Y is formed from Z, namely [ZI+o = [Zl-o = $ M - o [Ag+l+, = [&+I, - P I - 0 . [Ag] = - 1502.578 Phase Response of the Explodator However, TY] (i.e. y)? remains zero for a period Tuntil all the Ag+ has reacted. This time is reached when [Ag+l, = [Ag+I+,+[ [Ag+ldt = 0 (4) or [Ag+]+, = JOT 150zdt. (5) Since we had z = -502 clearly, the corresponding [Z] is [zl = [a+, exp (- 502) [&+I+, = 3[Zl+,( 1 - exp ( - 50731 (6) and from eqn (5) and (6) one obtains from which Tcan be calculated since [Ag+]+, and [Z]+, are known from eqn (2) and (3).Moreover, it can be seen from the last equation that the limiting case T+ 00 is associated with the relation By substituting eqn (2) and (3) in eqn (7), one obtains the critical Ag+ perturbation as at or above which the excess Ag+ cannot be consumed within a finite time. Hence, if one adds the critical or a larger amount of Ag+ ion and k = 0, the original system of differential equations (D') will never be applicable after the perturbation. If the pertuybation is smaller than critical, the dynamics of the system are described by the equations (DP) for a period of T, until the excess Ag+ is consumed. Thereafter, the equations (D') are again applicable, and the phase shift of the next maximum in [XI can be computed.(3) [Ag+], 3 2[yl~,+3[Z]~, As we have just discussed, this is the case of critical or super-critical perturbation. Without the limiting reaction (L 1) the phase shift is +m since in the absence of Y the [XI to be monitored will autocatalytically and monotonously increase ad infiniturn. Discussion With the given parameters (a = 100, b = 1, F = 0.01, k = 0) and a perturbation' of [Ag+], = 0.1, the sum of the numerical values 2y+ 32 = 2ry]-, + 3[z]-, is smaller than the perturbation in the interval 6.2 < t < 8.4 of the phase t of the stimulation (or on the relative timescale, with the period of the unperturbed limit cycle taken as unit: 0.653 < t < 0.884). Thus, the perturbation is supercritical in this interval, where the initial phase shift should be infinitely large. If an even larger amount of Ag+ ([Ag+], = 10) is used, the perturbation is subcritical only in the interval of 0.6 < t < 1 (or on the relative scale: 0.063 < t < 0.105).Hence, infinite phase shift should be found outside this interval. Although our results are in disagreement with the previous simulations' that used a different numerical approach, nevertheless we reach the same conclusion : with this specific set of parameters the Explodator cannot reproduce the experimentally observed phase response. Now let us examine a different set of parameters, with k > 0 in eqn (D). This corresponds to the inclusion of reaction (L 1). The limiting reaction (L 1) prevents the occurrence of infinite phase shifts since it offers a route for the production of Y from X [A~+I+, = 3[z1+,.(7) [Ag+lp, wit = 3[Zl-, + 2ryl-0 -f To preserve the correspondence between the numerical values of real and dimensionless concentrations, the symbol [ ] is retained for dimensionless concentrations.M. Eszterle, Z . Noszticzius and Z . A . Schelly - I phase of perturbation I I I I I I i I I I 579 a - c 5 f c a // - / / / / - / / / / - / / - / / / 0.5 0 -0.5 Fig. 1. Simulated phase response of the Explodator with (L 1) and (L 2) limitations, a = 2, b = 0.2, F = 0.01, k = 0.001. [Ag+], = 6. Fig. 2. Simulated phase response of the limited Explodator with (L 1) and (L 2) limitations. With the parameters given in this figure the experimentally observed excitable branch* is reproduced.a = 0.2, b = 0.02, F = 0.001, k = 0.01, [Ag+], = 15 (-), [&+I, = 45 (---).580 Phase Response of the Explodator which ultimately consumes the excess Ag+, Naturally, if excessive Ag+ perturbation is used, again, a new system of differential equations different from (D) must be applied so long as the excess Ag+ is consumed. However, with k > 0 this will always ensue, and the finite initial phase shift of the next [XI maximum can be computed. Clearly, the shape of the simulated phase response diagram is a strong function of the numerical values of the parameters used. For example, with a = 2, b = 0.2, F = 0.01, k = 0.001 and [Ag+], = 6 the phase response is depicted in fig. 1. The experimentally observed excitable branch of the BZ reaction2 can be reproduced by the limited Explodator, for example, with the following set of parameters: a = 0.2, b = 0.02, F = 0.001, k = 0.01 and [Ag+], = 15 or [Ag+], = 45 (fig.2). Since the actual rate constants of some of the reactions involved in chemical oscillations are usually not known, only estimated or trial values of the parameters can be used in the simulations. Hence, phase response relationships should be used with caution for diagnostic purposes. An agreement between experimental and simulated responses is necessary but not sufficient condition for the acceptability of a model considered. Recently, for instance, Dolnik et al." found the Brusselator12 superior to the Oregonator in predicting the phase advance observed experimentally in the BZ reaction following bromide ion perturbation, even though the Brusselator had not been designed to model specifically the BZ reaction. Consequently, the phase-response method can be viewed as a useful but not a crucial diagnostic test. This work was partially sponsored by the R. A. Welch Foundation. Acknowledgement is made to the donors of the Petroleum Research Fund of the American Chemical Society for additional support. References 1 P. Ruoff, J. Chem. Phys., 1985, 83, 2000. 2 P. Ruoff, J. Phys. Chem., 1984, 88, 2851. 3 A. M. Zhabotinskii, Dokl. Akad. Nauk USSR, 1967, 157, 392. 4 Biological and Biochemical Oscillators, ed. B. Chance, E. K. Pye, A. K. Gosh and B. Hess (Academic 5 R. J. Field and R. M. Noyes, J. Chem. Phys., 1974,60, 1877. 6 Z. Noszticzious, H. Farkas and Z. A. Schelly, J. Chem. Phys., 1984, 80,6062. 7 R. J. Field and R. M. Noyes, Faraday Symp. Chem. SOC., 1974, 9, 21. 8 P. Ruoff, Chem. Phys. Lett., 1982, 90, 76; 1982,92, 239; 1983, 96, 374. 9 Z. Noszticzius, J. Am. Chem. SOC., 1979, 101, 3660. Press, New York, 1973). 10 Z. Noszticzius, H. Farkas and Z. A. Schelly, React. Kinet. Catal. Lett., 1984, 25, 305. 11 M. Dolnik, I. Schreiber and M. Marek, Physica, 1986, 21D, 78. 12 I. Prigogine and R. Lefever, J. Chem. Phys., 1968, 48, 1695. Paper 71630; Received 9th April, 1987