Analysis of a mechanism of the chlorite-iodide reaction - The Journal

Mar 1, 1993 - Ian Millett, William Vance, and John Ross. The Journal of ... István Lengyel, Jing Li, Kenneth Kustin, and Irving R. Epstein. Journal o...
0 downloads 0 Views 3MB Size
Download PDF
J. Phys. Chem. 1993.97, 2851-2862

2851

Analysis of a Mechanism of the Chlorite-Iodide Reaction7 Peter Strasser, Janet D.Stemwedel, and John Ross’ Department of Chemistry, Stanford University, Stanford, California 94305 Received: September 30, 1992; In Final Form: January 5, 1993

We present a detailed analysis of the most recent mechanism of the oscillatory chlorite-iodide reaction proposed by Citri and Epstein. First, we use stoichiometric network analysis (SNA) to determine the complete set of major reaction pathways (extreme currents) of the mechanism under CSTR conditions. Two qualitatively different types of unstable extreme currents are predicted within the mechanism. At three different steady states the dominant extreme currents (EC) are singled out. Subgroups of elementary reactions (“dynamic elements”), dominant at different times during a cycle of oscillation, are identified. Suggestions for simplifications of the mechanism are given and a simplified oscillatory mechanism is extracted which accounts for the oscillatory behavior of the full mechanism. The SNA approach constitutes an operational prescription for the analysis of chemical reaction mechanisms, since the major steps of the analysis are based on computational routines. Second, the bifurcation behavior is investigated with particular attention to transitions from nonoscillatory to oscillatory regimes. We calculate one- and two-dimensional bifurcations with axes of experimentally convenient constraints. The roles of the chemical species in the mechanism are predicted based on the numerical results of the SNA analysis. We calculate the phase response of the mechanism as obtained from pulse perturbations of HOC1 and C102-. The calculations are supplemented by descriptions of feasible experiments to check the predictions. In the succeeding article bifurcations are determined experimentally and compared with the predictions presented here.

I. Introduction Although the complexity of the chlorite-iodide reaction has been known almost 30 years,’ the reaction was recognized about 12 years ago as the simplest oscillator of a whole family of oscillatingsystems.2 It has been used in connection with a variety of experimental investigations, which include studies on propagating waves in excitable media,3 stirring and mixing effects under CSTR conditions,4~5sustained pattern formation: and critical slowing down.’ In 1985 Epstein and Kustins set up a reaction mechanism for the chlorite-iodide oscillator, reported across-shaped phasediagram (XPD) in the ([Cl02-]o,[I-]o)-plane, and focused on the agreement of the experimental and calculated variations of the concentrations of measurable species with time. Two years later Citri and EpsteingJorevised the first suggestion for the mechanism by reducing the original model considerably. The batch behavior and the hysteresis behavior of the new mechanism were found to be in considerably better agreement with experimental results than those of the earlier mechanism. Olsen and Epstein” recently identified dimension two points within the parameter space of the Citri-Epstein mechanism in an article emphasizing the advantages of a bifurcation analysis over numerical integration of an initial value problem. Our study is based on the revised proposed mechanism for the chlorit&xlide reaction set up by Citri and Epstein? including the most recent adapted ratecoefficients?J2 Section 111is focused on an analysis of the dynamics of the model and for that we use stoichiometric network analysis (SNA).I3 According to this approach, we investigate the dynamical system by looking for its fundamental steady-state reaction routes called extreme currents (see Appendix). The complete set of extreme currents (EC) and their stabilities are calculated for the chosen reaction mechanism. Two qualitatively different types of unstable extreme currents are predicted to be present within the mechanism. The dominant extreme currents for three different steady states are singled out by means of a decompositionprocedure. These extreme currents constitute simplified submechanisms that provide insight into

’This and the following manuscript are dedicated to our colleague Harden McConnell in recognition of his great contributions and devotion to science. OO22-3654/93/2097-2851$04.OO/0

the origin of stability of the steady states. They enable us to extract a simple oscillator from which we identify subgroups of elementary reaction steps that determined the predicted dynamical behavior. In particular we identify the dominant subgroups of elementary steps (also called “dynamic element~”’~-I~) that are responsible for the limit cycle oscillationsof the mechanism near the Hopf bifurcation. With these identified dynamic elements we discuss the oscillations within the Citri-Epstein mechanism. Furthermore the extracted simple oscillator provides useful suggestions for possible simplifications of the full model. The computational method of finding the complete set of extreme currents and singling out the dominant ones constitutes an operational prescription for an analysis of chemical reaction mechanismsthat can easily be automated and transferred to other chemical systems. In addition to the predictions of the SNA approach, this study provides in section IV predictions of further numerical calculations,emphasizingglobaldynamical featuresand species-specific behavior, which can be checked by feasible experiments.I6 In section IV.l we extend the bifurcation studies undertaken in ref 11by using experimentally accessible bifurcation parameters and by investigatingthe nature of the transition between nonoscillatory and oscillatory regimes more closely. In the succeeding article bifurcations for this reaction are determined experimentally and compared with the predictions presented here. In section IV.2, based on the results of ref 21, we predict the roles of the chemical species in both types of unstable EC and indicate feasible experimentalprocedures to check the predictions. In section IV.3 we investigatethe phase response behavior of the oscillatoryregime by applying single-pulse perturbations of different strengths, and we present calculated phase transition curves together with suggested experiments. 11. Citri-Eptein Mechanism and Numerical Methods

Table I shows the set of elementary steps of the proposed reaction mechanism for the chlorite-iodide reaction. Table I1 gives the values and expressions of the rate coefficients used for all calculations. The dependence of the rate laws of reactions [ 11, [2], and [6] (Table I) on the concentration of H+is included Q 1993 American Chemical Society

2852 The Journal of Physical Chemistry, Vol. 97, No. 12, 1993 TABLE I: Elementary S t e p of the Citri-Epsteio Reaction MecMm"

H+ + CI(II1)

-

+ I-

HOCl

+ HOI

H++HOI+I---I,+H,o

+ H,O+ I-+ HOI + H+ HCIO, + HOI HOCl + HIO, H O C l + IHOI + C1H+ + HIO, + I- 2 H O I 2HOI HIO, + I- + H+ 2HI0, HOI + IO; + H+ HIO, + HOI I- + IO; + 2H+ HOCl + HIO, C1- + IO; + 2H+ I,

-

-

-

-

-

-

-

(1)

(2)

(3) (4) (5)

(6) (7)

(8)

(9) (10)

'The two reversible reactions ([2] and [6]) are written as two irreversible elementary steps ([2], [3] and [6], [7]).

TABLE II: Values and Expressions of tbe Rate Coefficients of the Cia-Epstein Reaction Mechanism k, = ( ( c , + c,[H+l)[H+l) [H+l + c3 c,

= 4.6 M-ls-'

c,

= 31 M-, s-'

c3 = 3.3 X

M

C4[H+1 )

k2 = (( [H'] + 0.002 M) c4 = 5 x 109 M-l S-I c4c5

k3 =

(([H'J + c,

0.002

M)

= 3.13 X lo-" M2

= 5.1 x 10-3 s - ~

k, = 1.4 X 10' M-' s-' k, = c,[H+] c,

= 1 x 106

k, = 25 M-l s-' k, = 3 X lo3 M-' s-' k, = 2.3 X 10, M-' s-' k,, = 1 X IO3 M-l s-' in the rate coefficients of Table 11. The reaction equations of Table I are always referred to by numbers in square brackets. The dynamical variables appearing in the differential equations are Cl(III), which constitutes a common variable for HC102and C102-, HOCl, H102, HOI, I-, and 12. The common variable for HC102 and ClOz- is justified by the fast equilibrium between both species. The pH is at 2.05.

Strasser et al. The numerical investigations of section 111, based on the concepts of stoichiometric network analysis (SNA)," are performed with Fortran routines17 designed for this purpose; these routines include the calculation of the complete set of extreme currents of a reaction mechanism (routine SNACONE), the identification and localization of unstable feedback cycles of a chosen reaction network (SNABETA), the calculation of the eigenvalues and of the concentration shift matrix (defined in section IV.2a) of a chosen CSTR network (SNAEVLS), and the decomposition of a chosen steady-state current into the closest basis of extreme currents which provides the dominant extreme currents (SNADCMP). All steady-state concentration values needed are calculated with the AUTO continuation package.', The Appendix gives a brief overview of the terminology and the basic ideas of SNA used in this study. The numerical bifurcation analysis presented in section IV.1 is performed by means of the AUTO continuation packagela with double precision used throughout. AUTO provides the continuation of branches of steady states and periodic solutions in one and two parameters. In the user supplied routine FUNC the time derivatives of all six dynamical variables are transformed into their logarithmicvalues before they are passed to AUTO this is necessary to avoid problems of convergence. The CSTR condition was modeled by adding the term

+kO([Xilo- [X,])[M/s]

i = 1, ...,6; M = mol/L

(1)

to the system of differential equations obtained from the elementary steps of the reaction mechanism. Thequantities [Xilo, [X,],and ko denote the input feed concentrations and the instantaneous concentrations of the species and the total flow rate,19 respectively. The initial stepsize for the units chosen is 0.01 and 5 X 1 6 ' for the one- and two-parameter continuation, respectively. We choose the total flow rate to be the first continuation parameter throughout. As the second parameter we choose [C102-]0 and the ratio [C102-]0/[1-]0. All tolerances are set to 1V. To obtain an initial guess of a steady state for the one-parameter continuation, we use the b o d e routine" for stiff systemsof ordinary differentialequations, choosing the initial flow rate ko to be 1.O s-I and integrating until the steady state is reached. In section IV.2 we use the sign pattern of the concentration shift matrix32provided by the cited Fortran routine that calculates the eigenvalues, in order to identify the role of the chemical species within chemical reaction networks. The Lsode routine for ordinary differential equations20 is also used to calculate the phase response curves in section IV.2.

III. SNA Approach for tbe Analysis of Steady-State D~MIII~CS under CSTR Conditions 1. Set of Extreme Currents of tbe CiM-Epstela Mecbmism. Our objective is to determine the complete set of fundamental reaction pathways of the reaction mechanism, called extreme currents (EC), and their stability. Eiswirth et al.21 first looked for extreme currents of the Citri-Epstein model close to the Hopf bifurcation by applying stoichiometric network analysis.13 We use a systematic computational way to obtain the full set of EC.I7 The procedure to determine the complete set of EC and their stability is implemented in the r6utines SNACONE and SNABETA, respectively. The terminology of SNA, used below, is described in the Appendix. The total number of EC of the Citri-Epstein reaction network is determined to be 39,of which 7 are unstable. The set of EC is obtained without using any rate coefficients. Since there are 18 reactions (10 chemical reactions plus 8 pseudoreactions) and no conservation conditions between the six variables, the EC form the edges of a 12-dimensional current cone. Figure la-c shows the network diagrams of the unstable EC that have the biggest percentage contributions, as we shall show when we discuss the decompositions of our three steady states. The EC in Figure la

The Chlorite-Iodide Reaction. 1

The Journal of Physical Chemistry, Vol. 97, No. 12, 1993 2853

a

d

(L

4L

16

9

t

1

4

1

Figure 1. Network diagrams of the most important extreme currents of the Citri-Epstein model mechanism. The species of the mechanism that are not dynamical variables (H+, Cl-, and IO3-)are not shown. Note that HC102 and C102- are equivalent. EC (a)-(c) are unstable and of type 1, (d) is unstable and of type 2. EC (e), (f), (g) are stable. The numbers next to the reactions refer to the elementary reactions in Table I. The number of feathers or barbs of a reaction arrow at a species correspond to the stoichiometric coefficients of the species in the respective reaction. The number of left feathers corresponds to the kinetic order of the species in the respective reaction.

is that reported by Eiswirth et al.21Figure Id shows an additional unstable EC of the model which represents the four remaining unstable EC since their topology is very similar to the EC in Figure Id. Figure le-g show the stable EC that have the biggest contributions to our steady states. In each current network diagram the number of feathers and barbs of a reaction arrow at a species represent the number of molecules of the species that are consumed and produced in the reaction, respectively. Since each EC fulfills the steady-state requirement, the total number of feathers at a species must equal the number of barbs. The number of left feathers corresponds to the kinetic order of the species in the respective reaction. The numbers next to the reactions refer to the numbering in Table I. EC la-c are of one type of topologically related unstable EC (referred to as type 1); similarly EC Id and the ones not shown but similar to Id are of another type, referred to as type 2. On combining either type of unstable EC with appropriate stable EC (e.g., to account for outflow pseudoreactions), choosing a total flow rate and looking at theeigenvaluesof the Jacobianof theresulting CSTR networks, we can show that either type can account for oscillatorybehavior. The procedure of evaluating the eigenvalues of an arbitrarily chosen network under CSTR conditions is implemented in the Fortran routine SNAEVLS.1' In addition to the eigenvalues,

SNAEVLS calculates the rate coefficients corresponding to the chosen network and chosen total flow rate. The values of the calculated rate coefficientsdiffer significantlyaccording to which type of unstable EC is used as the source of instability. Hence by means of SNA it is very simple to find different sets of rate coefficients, physically meaningful or not, for which a reaction network exhibits oscillations. As will be discussed further in section 1V.2, it is important to note that the chemical species play different roles2I in the two types of unstable EC. The instability of EC of type 1 is mainly based on the interplay of I-, HOI, and HI02, whereas that in the EC of type 2 is based on species I-, HOI, and ClOz-. HOI and HIOzform the autocatalytic cycle of the EC of type 1, and therefore they are called "autocatalytic species";zI I- and HOI have the same function in EC of type 2. The respective third variables, I- or C102-,play the role of the so-called "exit species",2' illustratedin Figure ld,since they consumeoneoftheautocatalytic species via an "exit reaction" producing inert products. It follows that two qualitatively different sources of oscillatory behavior are present within the Citri-Epstein mechanism, with the species playing different roles. Since we assume that species play only one role for a given set of rate coefficients, we expect the oscillations to be based on only one type of unstable EC. In

2854 The Journal of Physical Chemistry, Vol. 97, No. 12, 1993 EC2

Strasser et al.

TABLE IIk Decompositioa of Three Steady States of tbe Citri-Epstein Mechanism into Extreme Curreots~ EC (fig) % contribution stability (a) Stable Steady State on the High Iodide Branch 90.0 stable HOI, I equil 5 .o stable le 3 .O unstable la 1 .o stable If stable 1 .o 1g 3X10-4 unstable lb 1 x 10-7 unstable IC (b) Unstable Steady State Very Close to the Hopf Bifurcation Point on the Oscillatory Side HOI, I equil 97.0 stable la 2.0 unstable 1 .o unstable lb If 2 x 1v3 stable h3 3 X 1V3 stable 4x10“ stable le IC 2 x 10-6 unstable (c) Stable Steady State on the Low Iodide Branch 99.9 stable HOI, I equil 5 X 1V3 unstable la unstable lb 3 X 1V3 3 X 1V7 stable If 5 X 1V9 stable le 1 x 10-9 stable 1g 3 x 10-10 unstable IC ~~

Figure 2. The current Vdcmp is chosen to be decomposed into the extreme currents (edges) ECl, EC2, EC3, and EC4 of a three-dimensional closed steady-state cone. The closed steady-state cone shown is obtained from the original open cone by intersection with a plane of dimension 1 (here two-dimensional). The intersecting plane is indicated by the white basis of the upside down “cone pyramid”. Since the closed cone shown is not a simplicial cone (four edges instead of three) the decomposition is not unique, i.e., there is no unique vector j. To make the decomposition unique three out of the four EC have to be selected.

fact, when we discuss the decomposition of the unstable steady state close to the Hopf bifurcation, with the rate coefficients given according to Table 11, we shall see that only one of the two types of unstable EC is predominantly involved in the oscillations of the mechanism. The set of rate constants chosen alone determines which type of unstable EC is dominant. It follows that the roles of the species are dependent on the rate coefficients and therefore may vary within the parameter space. 2. DeterminationoftbeDominantECofSteadyStates. Having found the set of EC of the Citri-Epstein mechanism, we determine the dominant EC for each of three steady states (a stable steady state at a high value of ko on the high iodide branch, another stable steady state at a very low value of ko on the low iodide branch, and an unstable steady state at an intermediate value of ko and close to the Hopf bifurcation point on the oscillatory side) by decomposing the corrrespondingcurrents (uhigh, ulow,and uHopr) into the EC of the mechanism. The three u are obtained from the steady-statevalues of the six dynamical variables (calculated by means of a continuation routine’*)and from the values of the rate coefficientsaccording to eq A1 in the Appendix. The steady state current which is to be decomposed is referred to as udcmp. Figure 2 shows udcmp for a current cone e, of dimension 3 which has four edges, i.e., extreme currents. The following decomposition procedure of the current vector udcmp describes the determination of the dominant EC (implemented in the routine SNADCMP): If there are no conservation constraints of the variables present, as in the Citri-Epstein mechanism, the dimension of the current cone e,,containing all the steady state currents, equals (r - n), where r is the total number of elementary steps (including the pseudoreactions) and n is the number of dynamicalvariablesof the mechanism. Usually, however, the number of extreme currents (i.e., the number of one-dimensional edges of e,) is greater than (r - n), as shown in Figure 2, where (r - n) equals 3 and the number of extreme currents equals 4. Therefore the vector udcmp can not be uniquely expressed as a linear combination of the EC. According to eq A7 this is equivalent to the fact that there is no unique weighting vectorj that would correspond to the current udcmp. To make the decomposition unique, (r - n) extreme currents of the full set of EC must be selected, and these must be those which are closest to the current udcmp, since they are topologically most similar to Udcmp. The distance between udcmp and the EC is determined as follows: e, is intersected with the hyperplane &j = 1, which is of dimension one. The (r- n) closest EC are selected according

* The first column refers to the network diagrams of Figure 1, the second column gives the percentagecontribution, the third column shows the stability of the EC in column 1. to the distance of the intersection points of the EC vectors to the intersection point of u d m p . Simultantously two additional conditions must be fulfilled by the selected set of EC: first, all the EC must be linearly independent and thereby span a simplex; second, the simplex must contain the vector u d m p in its interior. The selected EC can be considered as the edges of a new cone e’,, again of dimension (r - n), within the old cone e,. The linearly independent EC form a basis of e’,, which is also referred to as a simplicial cone,13 since it has a minimal number of edges but the same dimension as the original cone. The hyperplane &ut = 1 splits the simplicial cone e’, into a closed portion containing the origin (see Figure 2) and an open portion. In two dimensions the closed portion of the simplicial cone has the shape of a triangle, in threedimensions,a tetrahedron. Thedecompositionof a current u d m p which lies in the interior of the simplicial cone e’, into the edges (extreme currents) is now uniquely determined by nonnegative weighting coefficients;there is a unique vectorj according to eq A7 with nonnegative components. For different U d m p the simplicial cone e’, must be recalculated, since the closest (r- n) extreme currents fulfilling the required conditions may change. By means of the simplicial cone we find a unique vector j for each udcmp. The componentsji of the vector j represent the relative contributions of the selected extremecurrents to the current Ddmp. The jj are converted into percentage contributions. The EC correspondingto the biggestjiis mainly responsiblefor the stability of the current u d m p and is therefore the dominant EC. Table 1IIa-c shows the decomposition of the steady state on the high iodide branch, near the Hopf bifurcation point, and on the low iodide branch, respectively. The entries of the first column refer to the network diagrams of Figure 1. The percentage contributions of the EC in the second column correspond to the respectivemagnitude of theji. The third column gives thestability of the EC. We see that all three steady states are dominated by the 12, I- stable equilibrium EC. Equilibrium EC, however, are not of interest for far-from-equilibrium phenomena like the high iodide branch or the Hopf bifurcation. Therefore let us consider the EC le which is the most dominant, nonequilibrium, stable extreme network of Table IIIa (high iodide state). This extreme networkconsistsofthechemical reactions [l], [2],and [5] (Table I). This result is consistent with results in ref 9, where the limit

The Chloritc-Iodide Reaction. 1 cycles of the Citri-Epstein mechanism are interpreted in terms of elementarysteps by analysisof phase portraits: In ref 9 reactions [l], [2], and [5] are found to be dominant for the slow nonautocatalyticconsumptionof I-. Together with the high inflow rate of I- (EC If) on the high iodide (flow) branch two opposing forces are present which stabilize the steady state at a high concentration value of I-. Despite its contribution of 3% the unstable EC l a is not able to cause an overall instability since the high iodide branch is stable. There are no clear rules as to what kind of unstable EC must participate and to what extent they must contribute in order to cause instability. Stable EC lg, which is based on the CSTR arrangement of the system like EC 1f, enhances thestability of the network further. Thecontribution of the next unstable EC 1b is negligible. The unstable steady state near the Hopf bifurcation shows a very different contribution pattern as can be seen in Table IIIb: The trivial equilibrium EC is increasingly dominant. Now the unstable EC l a becomes very important and is, together with the unstable EC 1b, responsible for the instability of the steady state. All other stable EC contribute very little to the overall network. It follows that given the rate coefficients of Table I1 the EC of type 1 are responsible for the oscillations of the model close to the Hopf bifurcation, and this result confirms the correctness of the result in ref 21, EC la consists of the chemical reactions [2], [4], [5],and (61 andcontainsanautocatalyticcyclein HOI. This is again consistent with results of ref 9 where reactions [4], [5], and [a] were found to be crucial for the autocatalysis of HOI which is necessary to bring the system back to a high (low) HOI (I-) concentration. Finally let us consider the decompositioninto dominant extreme currents of the stable steady state on the low iodide branch, Table IIIc. The low iodide branch is the bridge from the far-from-equilibrium behavior to the equilibrium behavior, since it is the only stable attractor when the CSTR arrangement approachesflow ratesvery close to zero. Accordingly we expect the steady state to be determined by equilibrium rather than far-from-equilibrium dynamics. Table IIIc confirms that the contribution of the equilibrium EC is almost 100%. The unstable and stable EC lose importance by 3 orders of magnitude. 3. Simplified Oscillatory Mechanism and Its D~MIO~CS. In this section we first focus on the determination of a simplified oscillatory network. Since the oscillatory network is based on stable and unstable EC of the full mechanism that are dominant near the Hopf bifurcation point (thereby all unstable EC participatingin thesimplifiedoscillatorareoftypc l), thenetwork can be viewed as an extracted oscillator which forms the basis of the oscillations of the full mechanism. Within this network we then identify subgroups of elementary reactions, which dominate the oscillations at different times during a cycle of oscillation. These. subgroups determine the dynamics of the full Citri-Epstein mechanism near the Hopf bifurcation. As reported by Aguda and Clarke14J for the peroxidawxidase model, dominant subgroups of elementary reactions, which are called “dynamicelements” in refs 14and 15, can be determined directly from the full mechanism of the oscillator. However, the dominating EC, which constitute simplified subnetworks, can ease the determination of dynamic elements considerably if the full mechanism is very complex. We show how the determination of ‘dynamic elements” of the Citri-Epstein mechanism is performed using the simplified submechanisms(EC) rather than the full model. Our objectiveis to find the simplest combinationof thedominant EC (Le., the simplest network) that exhibits oscillations under CSTR conditions, and we use the procedure implemented in the routine SNAEVLS:” The value of the total flow rate ko for the CSTR networkischosentobeequal to that at the Hopf bifurcation point of the full mechanism, thereby preservinga Hopf bifurcation point of the simplified mechanism near the chosen value of ko. In addition, a vector ] is chosen that determines which EC contribute to our network and to what extent. We choose the vector] that has been provided by the routine SNADCMP after

The Journal of Physical Chemistry, Vol. 97, No. 12, 1993 2855 the decomposition of the unstable steady state near the Hopf bifurcation point. According to the 12-dimensionalsteady-state cone, the vector j contains 12 nonzero componentsj i (see section 111.2). Since the vector j corresponds to a steady-state reaction rate vector v(see eq A7). and since CSTR conditions (Le., the same total outflow rate for each species) are required, the steadystate values of all six speciescan be uniquely determined by solving eq A1 of the six outflow reactions for the six steady-state concentrations Xio(it is assumed that the first 12 components of the reaction rate vector v, u1-012, correspond to the 10 chemical and the two inflow reactions, whereas the last six components of v, ~ 1 through 3 ~ 1 8 ,correspond to the six outflow reactions):

Xi, = ui+,2/ko Xi, = (E,)i+,z/ko

i = 1, ..., 6 i = 1, ..., 6

The columns of the matrix E are the extreme currents of the original network. The rate coefficients kr of the 10 chemical reactions and of two inflow reactions can now be determined by solving eq Al:

k, = q/nq;

1 = 1, ..., 12

(3)

i= 1

where ~~l is an element of the kinetic matrix and denotes the kinetic exponent of the ith species in the Ith reaction. The rate coefficients are provided by SNAEVLS given ko and a vector j. Step by step the initial vector j is simplified by zeroing nonzero components j i , thereby decreasing the number of participating EC. The oscillator consists of fewer and fewer reactions, making the network simpler and simpler. Since we want the simplified oscillator to be based on the most dominant EC, the smallest nonzeroj i are the first to be set to zero. Furthermore, thej, that are set to zero must be chosen such that the CSTR conditions (i.e., the presence of all six outflow reactions) are fulfilled for the simpler mechanism. SNAEVLS indicates whether the chosen j still fulfills the CSTR-conditions. After each zeroing of a j i the eigenvalues of the remaining network corresponding to the simplified vector j are monitored using SNAEVLS. A simplification is considered to be successful if the simplified mechanism corresponding to the simplified vector ] still exhibits stable oscillations, i.e., positive real part and nonzero imaginary part of two eigenvalues, and negative real part of all other eigenvalues. The procedure of setting j i equal to zero is repeated until no j i can be set to zero without making the oscillation cease. Following this procedureonly thechemical reactions [2], [3], [4], [5], [a], and [8] (plus all eight flow reactions) are found to be necessary for a CSTR-network to exhibit stable oscillations at the chosen value of ko. The rate coefficients of this network do not differ appreciablyfrom thoseinTable 11,sincethevaluesof thedominant j i remain unchanged. It follows that a CSTR mechanism containing these six chemical reactions constitutes a simplified oscillator (henceforth referred to as ‘oscillator 1”) when quantitative agreement of the new and the original set of rate coefficients is required. Oscillator 1 is built on EC la,b,c,f,g, the equilibrium EC consisting of reactions [2] and [3], and two other stable EC; out of the 12 nonzero components j i of the original vector] there are eight componentsleft in vector jorl corresponding to oscillator 1. The role of the essential species within oscillator 1 corresponds to that within the unstable EC of type 1 since all the unstable EC contained in the linear combination are of type 1 (EC la,b,c). We have extracted a simplified mechanism from the full mechanism, which exhibits essentially the same local bifurcation behavior for the same values of the rate coefficients, with the species playing the same role. Since, however, we are only interested in a qualitative analysis of the oscillations (Le., we want to extract the elementary steps of the full mechanism that are necessary for stable oscillations), we also want to eliminate

2856 The Journal of Physical Chemistry, Vol. 97, No. I2, I993

those elementary reactions of oscillator 1 that do not affect the oscillations qualitatively but instead are merely responsible for thequantitative agreement between thesimplified and the original set of rate coefficients. To this end, we now allow a change of the rate coefficients during the search for a further simplified oscillatory mechanism. The condition of a Hopf bifurcation at the value of the total flow rate chosen, however, is still in place. Furthermore, we require that the species play the same role in the simplified mechanism. Thereby we achieve a clearer qualitative insight to the origin of the limit cycle oscillations at theexpenseof somequantitativedetails. We start with thevector jWl and vary the values of the j i smoothly rather than zeroing them abruptly. Thereby we try to zero additional components j,. If j , of high values (corresponding to dominant EC) need to be changed in order for the simplified mechanism to meet the required oscillatoryfeatures, the newly calculated rate coefficients of thesimplified network may differ significantlyfrom the original ones. We simplify the vector jwI further such that as many reactions as possible are eliminated, but such that the remaining mechanism still exhibits a Hopf bifurcation near the total flow rate chosen. Stable oscillations and a Hopf bifurcation near the required value of the total flow rate can be obtained simply by superimposing EC la,b,f,g and two other stable nonequilibrium EC. The equilibrium EC consisting of reactions [2] and [3], and therefore reaction [3], is not needed. Only unstable EC of type 1 (EC 1a.b) are contained in the linear combination. Therefore, the role of the species corresponds to that of the full mechanism. The stable EC are necessary to stabilize the unstable EC, which simply lead to an explosion of concentrations when used alone. Furthermore, the stable EC ensure that the CSTR conditions are fulfilled,which means all six outflow reactions are present. Figure 3 shows the network diagram of the further simplified oscillator without reaction [3] (from now on referred to as 'oscillator 2"). The six outflow reactions are omitted. Oscillator 2 contains only the chemical reactions [2], [4], [5], [6], and [8], but according to the eigenvalues there are all sufficient qualitative features for stable oscillations present. All attempts to reduce oscillator 2 further, without making the oscillations cease, have failed. Note that Figure 3 is a steady state network, since it constitutes a linear combination of EC. However, since the coefficients of the linear combination are partly nonintegers, the resulting stoichiometries of oscillator 2 are not simple multiples of the reactions in Table I. Therefore, for simplicity the stoichiometry of each individual numbered reaction in Figure 3 is chosen to be the same as that in Table I. By looking at oscillator 2 (Figure 3) we are able to identify qualitatively two "dynamic elements" (DE), to which we refer as DE-1 and DE-2: DE-1 consists of two parallel autocatalytic cycles of HOI which are formed by reactions 141, [6] and [4], [ 5 ] , [8]and by species HOI, HI02 and HOI, HOC1, respectively. DE-2 consists of the inflow reactions of both I- and

c102-. The remaining reaction [2] is a very fast "exit-reaction", since it takes care of the fast depletion of either HOI or I-. Reaction [2] is reminiscent of what is called a 'switch" reaction in ref 15. Switcheshave been found to be of great importance in oscillatory systems.I4J5 Their characteristic phase portrait pattern is the extinction of one competing variable and the simultaneous explosion of the other provided that their rates of productions are not equal (see Figure 5 of ref 15). The phase portrait in the HOI, I-planeof concentration spaceshowsthis typical hyperbolicshape. Now we show how to use these identified dynamic elements to describe the oscillations. We find the following interplay of the dynamic elements responsible for the oscillations: Suppose the I- and the C102- concentrations are low whereas the HOI, HIO2, and HOCl concentrationsare high as predicted at a certain phase of the limit cycle by calculations of phase portraits. The inflow reactions of DE-2 increase [I-] and [CIOz-] steadily. Simultaneously reaction [2] (the switch element) lowers [HOI]. HOCl and HI02 are transformed into H 0 1 via reactions [ 5 ] and [a] and therefore are also depleted. The rate of depletion of HOI

Strasser et al.

I

I2

I DE-2

"t 1

E HI02

I

I

DE-1

Figure 3. The network diagram of oscillator 2 which was constructed by linearly combining EC 3a,b,f,g and two other stable EC. This oscillator is used to illustrate qualitatively the origin of the oscillations of the full Citri-Epstein mechanism. All the unstable EC contained in the network are of type 1, which has been found to be the dominant type near the Hopf bifurcation point. Thequilibrium EC, which is the only EC that contains reaction [3], need not be included in the linear combination in order to qualitatively account for oscillations. The set of ratecoefficients. however, for which oscillator 2 shows oscillations, differs from the original set given in Table 11. On including reaction (31 in the linear combination, a quantitative agreement of the simplified and the original set of rate coefficients is possible. It follows that reaction [3] affects the dynamics of the mechanism only quantitatively, not qualitatively. The six outflow reactions of the CSTR network are not shown. The numbers next to the reactions refer to Table I. Figure 3 constitutes a steady-state network for the purpose of clarity the stoichiometries of the reactions (number of feathers and barbs) are chosen to be as given in Table 11, and hence the number of feathers at a specie-s does not necessarily q u a l the number of barbs. The netwock can be considered as a submechanism extracted from the full mechanism, which is dominant near the Hopf bifurcation.

is governed by the inflow rate of I-. Extinction of HOI, however, is never reached completely sinceas soon as C102-has accumulated enough the autocatalysis of DE-1 is turned on and the rate of production of HOI is suddenly higher than that of I-. [HOI] has become very low at this point. Due to the very high rate of production of HOI at this stage, the switch element now causes a very fast drop in [I-] according to its phase portrait (Figure 5 in ref 15). I-, however, is crucial for maintainingthe autocatalysis via reactions [5] and [6]. Hence as soon as I- drops below a certain value, the autocatalysis is turned off and the rate of production of HOI drops considerably,whereas I- is still produced at a constant rate by the inflow reaction. Now the system is back in the initial state. It follows that DE- 1 and DE-2 are responsible for the increase in [HOI] and [I-], respectively. Both dominate the full mechanism at different time intervals. The depletion of either species is mainly taken care of by reaction [2]. On the basis of the analysis given above, we expect an increase of the period of the limit cycles if the ratio [ClO2-]o/[1-]o of the two input feed concentrations is lowered, since for smaller ratios it takes more time for CIOz- to build up the concentration that is necessary to turn on the autocatalysis of HOI. For increasing values of the ratio, however, we expect the following behavior according to our analysis: I- has less time to accumulate since the autocatalysis is turned on earlier; hence the maximal concentration value of I- becomes smaller. Simultaneously the minimal concentration value of HOI increases, since there is less time for it to be depleted. The result is a steady decrease of the amplitude of the periodic orbits, with increasing values of [ClOz-]0/[1-]~. The predicted behavior of the model, which is based on our qualitative explanationof the limit cycles, for higher

The Chlorite-Iodide Reaction. 1 and lower values of the ratio of the two input feed concentrations will be confirmed by calculations in section IV. 1. 4. Suggestions on Possible Simplifications of the Mechanism. Oscillator 1, which constitutes a simplified oscillatorymechanism extracted from the full Citri-Epstein mechanism, suggests reactions which may be neglected without changing the overall dynamics significantly.22Oscillator 1 has been found to exhibit very similar local bifurcation behavior for essentially the same values of the rate coefficients. Since reactions [l], [7], [9], and [lo] are missing in oscillator 1, we can conclude that they are of minor importance for the oscillations near the Hopf bifurcation point, which occurs at a medium value of the total flow rate ko. In section 111.2, however, reaction [ 11 was found to be of great importance for the nonautocatalytic consumption of high concentrations of I- and Cl02- at high values of k ~ .Therefore we conclude that for medium values of ko reaction [l] is only responsible for the transient behavior of the mechanism; the reaction simply decreases high initial concentrations of the flow species such that oscillations can occur. For high values of ko reaction [ 11regulates the highconcentrationsof theidlow species, and thereby makes the existence of a high iodide (flow) steady state possible. Calculations of the global bifurcation behavior of the mechanism without reaction [ 11 confirm our conclusions. As indicated in the previous section, reaction [3] is necessary for a good quantitative agreement of the simplified and the full mechanism but is qualitatively unnecessary for oscillations. Reactions [7], [9], and [lo] have been found to be of minor importance for small, medium, and high values of ko. Therefore we conclude that all three reactions can be neglected without changing the dynamics of the mechanism significantly. Calculations of the bifurcation behavior again confirm our conclusions. However, if all three reactions [7], [9], and [lo] are neglected, the overall stoichiometry of the chloriteiodide reaction cannot be reconstructed by the elementary steps of the simplified mechanism. At least one of the three reactions must be included in order to meet the required stoichiometry. From the above analysis the following possible simplifications of the Citri-Epstein mechanism are suggested: Two reactions of [7], [9], and [ 101can be neglected, thereby preserving the global and local dynamical features of the mechanism, the values of rate coefficients and the overall stoichiometry of the chlorite-iodide reaction.

IV. Predictions and Suggestions for Experiments

The Journal of Physical Chemistry, Vol. 97, No. 12, 1993 2857 16.5 (0

.-

U

15.5 c

0

14.5 C

1354

.

-3.3

.

.

-2.8

-2.3

-1.8

log [ k o l

Figure4. Bifurcation diagram in one dimension with koas the bifurcation parameter; the norm of species, eq 6, is plotted vs the logarithm of the total flow rate (inverse of residence time), ko. The values of [C102-]0 and [I-10 equal 3.187 X le3and 8.563 X le3M, respectively. Solid lines denote the stable high iodide (SSl)and the stable low iodide steady state (SS2), respectively. Thedashed linedenotes unstablesteady states. The dotted-dashedline (po) denotes the branch of periodic orbits. Further abbreviations are: snl, sn2 = saddle-node infinite period bifurcation point; snp = saddle-nodebifurcation point of periodic orbits; HP = Hopf bifurcation point.

The total flow rate ko is chosen to be the constraint for the one-parameter continuation. The following relations hold for ko and the input feed concentrations [ i ] of ~ species i:

D, i

ko=-=-

V

A.

[ilO= -Ci

D,

R

D/sl

(4)

f;:

= -Ci

VkO

[MI

(5)

i

with V, R , J , and Ci denoting the reactor volume, the residence time, the individual flow rates of C102-, I-, and H+ and their bottle concentrations, respectively. Figure 4 shows the one-dimensional bifurcation diagram on variation of the constraint of the total flow rate, which we briefly discuss in preparation for the two-dimensional bifurcation diagrams. The pH equals 2.05 throughout. The values of the input feed concentrations are 3.187 X le3and 8.563 X le3M for [C102-]oand [I-lO,respectively. We choose as they axis the norm of all six species concentrations, which is defined by

1. Bifurcation Structureof the Parameter Space. (a)Preuious 6 Work. Theoretical and experimental investigationsof the parameter space of the Belousov-Zhabotinski reaction were reported llull = [ ~ u ~ ] I / for 2 steady states by various a ~ t h o r s ~ more ~ - ~ ’than a decade ago. The predicted i= I cross-shaped phase diagram of the Briggs-Rauscher reaction was experimentally determined by Boissonade and DeKepper28by llull = !-$ku;(t) dr]’” for periodic solutions ( 6 ) varying two experimentally accessible constraints to obtain a twoi= I dimensional cut through parameter space. Similar experimental studies have been applied to the chlorite-iodide ~ y s t e m . s ~ ~ J ~ The ~ ~ ~two norms30defined by eq 6 are compatible in that they All the theoretical studies, however, are based on the original have the same value for Hopf bifurcations and saddle-nodeinfinite 13-step and only the input feed concentrations of chlorite period bifurcations in the bifurcation diagram, regardless of and iodide are used for the two-dimensional bifurcation plots. In whether the point is considered as a stationary or a timedependent 1991 Olsen and Epstein” studied the bifurcation behavior of the solution. The plot shows two branches of stable stationary states Citri-Epstein mechanism, identifying dimension two points SSl and SS2 which are of high and low iodide concentrations, within the parameter space of the mechanism. In the same year respectively. The stable branches are connected via two saddleRingland30 reported an interesting approach to test two BZ node infinite period bifurcations snl and sn2. The transitions mechanisms by predicting certain bifurcation features and between oscillatory and nonoscillatory regions, four main types comparing them with experimental results. of which were described by AndronoP and later conceptually applied to chemical ~ystems,25-2~J~~I are of particular interest. ( b )Numerical Results andSuggestions for Experiments. We The oscillatoryregion is indicated by the branch of periodic orbits calculate one- and two-dimensional bifurcation diagrams of the (po). Hysteresis between stable orbits and the stable branch of Citri-Epstein model9 with the convenient choice of constraints low iodide steady states (subcritical Hopf bifurcation) is predicted of the total flow rate ko,I9the input feed concentration [C102-]0 for the transition at low values of ko, whereas a saddle-node infmite and the ratio [C102-]0/ [I-10 and suggest experimental investiperiod bifurcation without hysteresis is predicted for the transition gations.

[

2858 The Journal of Physical Chemistry, Vol. 97, No. 12, 1993

1

/

t u

Y

-m

I -2.84 -4.0

Strasscr et al.

II

. .

/

ss 1

I I

. . . . :

-3.0

-2.0

-1.0

0.04

-3.0

-4.0

-2.0

-1.0

log [ k o l

Figure 5. Bifurcation diagram in two dimensions with koas the first and with [C102-]0as the second axis variable. The value of [I-lo equals that in Figure 4. SI, SS2 = region of high and low steady states; Osc = region of periodic orbits; Exc = region of excitability; snl, sn2 = curves of saddle-node infinite period bifurcations; hp = curve of Hopf bifurcations; snp = curve of saddle-node bifurcations of periodic orbits; swt = swallow tail (small area of tristability); C = cusp point; TB = Takens-Bogdanov point. The horizontal line symbolizes the location of Figure 4 within Figure 5: consider Figure 5 to be extended in a third dimension by a artesian axis I denoting the norm according to q 6. Figure 4 is obtained by intersecting the three-dimensional extension of Figure 5 with a twodimensional plane, which is parallel to the ( 2 , ko) plane, and viewing the intersection in direction of the [C102-]0 axis.

from oscillation to the branch of high-iodide steady states. Excitability is predicted for the system located on the branch of high iodide steady states for values of ko between those of the two saddle-node infinite bifurcation points snl and sn2 (see Figure 4). We extend the one dimensional bifurcation set of Figure 4 by choosing as a second constraint to be varied the input feed concentrationof CIOz-, leaving theinput feedconcentrationvalues of I- and H+ fixed. Figure 5 shows the two-dimensional bifurcation set in the (ko,CIOz-) plane. The regions of distinct dynamical behavior like oscillations (Osc), excitability (Exc), and stable steady state (SSl and SS2) are marked. Hp denotes the branch of Hopf bifurcation points, snl and sn2 the branch of steady-state saddle-node infinite period bifurcations, and snp the branchof saddle-nodesofperiodicsolutions. Other interesting features of the model are the change from s u b to supercritical Hopf bifurcation points for high and low values of [Cloz-]~and the swallow tail of the branch of steady-state saddle-node infinite period bifurcations, which usually is not obtained for a minimal oscillatory model. A number of predictions of the calculated bifurcation diagram can be tested by feasible experiments: For the required fixed values of the input feed concentrations of Iand H+, increase the individual inflow ratesf) by the same ratio. According to formulas 4 and 5 ko is increased whereas all the input feed concentrationsare left unchanged. This is a convenient way to monitor the experimental behavior parallel to the ko axis. To move up and down towards higher and lower values of [C1O2-]0, change the bottle concentration of C102-. An alternate way to monitor the bifurcation diagram in the vertical direction without changing bottle concentrations is to vary the ratio of [ClOZ-]oand [I-]o by changing the respective fr, leaving their sum constant. This corresponds to a vertical movement in the bifurcation diagram, with the movement along the abscissa as before. We calculate an additional bifurcation diagram (Figure6) using [ClOz-]o/[I-]~asthesecondcontinuation parameter. The regions of distinct dynamical behavior are denoted as in Figure 5. Since [I-IOis no longer constant, Figures 5 and 6 are not identical, although the two cuts through the parameter space are close to each other. The period and the amplitude of the oscillations are to be monitored by measuring the instantaneous I- concentration. Interesting dynamical phenomena are predicted by the model mechanism for 0.25 < [C102-]0/[1-]0 < 0.65,32asshown in Figure 6. The ratio [C~OZ-]O/[I-]O is denoted as y. We start at point (1) within the oscillatory region and approach the curve snl of

18i

8i 0

0 0

b

0.02

0.04

0.06

0.08

0.1

0.12

Distance from Bifurcation Point P O ~ - ] 0 ~ 1 ~ - 1W0 ~ O ~ ~ l

~ ~ [ ~ ~ l o ~ ~ i ~ Figure7. Oscillatory period as a function of distance from the bifurcation point for saddle-node infinite period bifurcation. Data calculated from the Citri-Epstein model at ko of 2.667 X lO-' s-I, plotted on axes of period versus difference between [CIO~-]o/[I-]ofor the oscillating point and for the bifurcation point.

saddle-node infinite period bifurcations along the horizontal line by increasing the total flow rate (e.g., along the line shown in Figure 6). The period increases, regardless of the value of y. Moving towards higher values of ko,a smooth transition to the high iodide state occurs at the intersection point (2) of the curve snl of saddle-node infinite period bifurcations and the horizontal line that corresponds to the c h a m value of y. The amplitude at the transition point is non-zero and the period becomes infinite. Figure 7 shows the exponential increaseof the period of oscillation upon nearing the bifurcation points4* There is no hysteresis in theneighborhood of point (2); on lowering &O the transition back to the limit cycle occurs again at the value of ko of point (2). On increasing ko further, the calculations predict a region of excitability that starts at the ko value of point (2) and ends for ko of point (3). The region of excitability decreases in width with decreasing values of y . For higher flow rates (point (4)) the stable steady state becomes nonexcitable. We start again at point (1) in Figure 6 and decrease ko toward the left border of the bifurcation diagram; a "hard transition" (subcritical Hopf bifurcation) to the low iodide state SS2occurs for 0.3 < y < 0.52. In Figure 6 the transition point ( 5 ) occurs at the intersection of the curve snp and the horizontal line. The hard transition is characterized by a nonzero amplitude and finite period at the transition point. At thesubcritical Hopf bifurcation of Figure 6 bistability between SS2and the limit cycle is predicted for the parameter region between the curve hpof Hopf bifurcations and the curve snp of saddle-nodes of periodic orbits; hence, on

The Chlorite-Iodide Reaction. 1 increasing ko,the transition back to the limit cycle occurs at point ( 6 ) instead of point ( 5 ) . This transition is again hard; amplitude of the oscillationsincreases abruptly at the bifurcation point ( 5 ) from zero to a nonzero value. For very low values of ko (point (7)) the model predicts a nonexcitable stable steady state. For y > y l and y < yz the transitions between SS2 and the limit cycle occur via a soft transition (supercritical Hopf bifurcation) when crossing the curve hp: The amplitude of the oscillations goes smoothly tozero, whereas the period remains finiteat the transition point. The reverse is observed for transitions from SS2 to the oscillatory region. There is neither excitability nor hysteresis near the supercritical Hopf bifurcation. The same type of transition is predicted for the upper border of the oscillatory region, reached by holding the total flow rate constant and increasing y . Note that, while excitabilityhas been reported in other systems before subcritica143 or ~upercritical~~ Hopf bifurcations, our calculations do not predict these phenomena. Moreover, since regions of hysteresis may be extremely the experimental distinction between a saddle-loop bifurcation and a saddle-node infinite period bifurcation is tentative at best. A number of features of the bifurcation diagram discussed here are confirmed experimentally in the succeeding paper; the region of excitability is substantially larger than that predicted.16 2. Identification of the Role of the Model Species. ( a ) Classification of the Role of Species.z1 On the analysis of 25 systems of abstract models and realistic chemical mechanisms Eiswirth et ale2'distinguish four categories of oscillators. The assignment of a system to one of these categories uses the role of species in the dominant unstable extreme current as a major criterion. The primary classificationof species is their distinction into essential and nonessential for a Hopf bifurcation. Three different kinds of nonessential speciesare defined. The remaining essential species are defined according to their role in the mechanism. Depending on the category of the oscillator, up to four classesof essential species are defined: autocatalytic species denoted as species of type X, "exit species" of type Y,"feedback species" of type Z and "recovery species" of type W. Autocatalytic species are characterized by their participationin the autocatalysis; exit species react with an autocatalytic species; feedback species are consumed by the autocatalytic cycle; and finally recovery species are usually produced by reactions of species of type Y and X, and thereafter they react again with Y species. According to ref 21 an operational test for a distinction of essential species is a constant increase or decrease of their inflow concentrations in a CSTR. A constant addition of a species shifts the steady state concentrations of all species in the mixture either to higher (+) or lower (-) values. Only the effect of the addition of each species on the steady-state concentration of the same species need be considered, Le., the diagonal element of the concentration shift matrix. An element mi&of this matrix denotes the steady-state concentration shift of the ith species caused by a small constant addition of the kth one. To make the assignment of the essential species unique, ref 21 indicates a second role-specific criterion: If a constant addition of a species is applied at a stable focus close to a Hopf bifurcation, the system can either begin to oscillate or remain stable. In the first case the system is destabilized (d), in the latter case stabilized (s) by the addition. Accordingly, if we apply a constant perturbation starting at the oscillatory side of a Hopf bifurcation, the oscillationsmay either cease (s) or remain (d). On the basis of calculations of the concentration shifts and the effect of perturbations on the stability ref 21 provides the followingresult (seealso Table IV, columns 2 and 3): autocatalytic species X and recovery species W are s+, the exit species Y are d+, whereas feedback species Z are s-. ( b )Numerical Results and Suggestionfor Experiments. For the purpose of identification of the roles of the species in the model we use a method33 based on the sign pattern of the elements of the concentration shift matrix including both diagonal and off-diagonal elements. The model predicts that the species of the

The Journal of Physical Chemistry, Vol. 97, No. 12, I993 2859 TABLE I V Essential Species of the Citri-Epstein Mechanism Shown in Column 1 As Predicted by the SNA Analysis, Their Types Shown in Column 2 As Predicted by Calculations of the Concentration Shift Matrix,3* and Their Experimental Behavior Symbolized in Column 3 As Predicted by Means of the Methods of Ref 21' essential species CIO*IHOI HI02 HOCl

predicted tY Pe

predicted effect

2

S-

Y

d+

X X X

S+ S+ S+

X denotes "autocatalytic species", Y denotes 'exit species", Zdenotes "feedback species". The symbols d (s) indicate a destabilizing (stabilizing) effect of the species on the Hopf bifurcation. The signs + and - are those of the diagonal element of the concentration shift matrix of the respective species. (I

unstable EC of type 1 (see Figure l a x ) play the following roles: Cl02- is a type Z feedback species, I- is a type Y exit species, whereas HOCl, HI02, and HOCl are all type X autocatalytic species. I2 is an inert product and therefore nonessential. The roles of the species in the unstable EC of type 2 (see Figure Id) are found to be as following: HI02 is a species of type Z, CIOzis of type Y,whereas I- and H 0 1 are of type X. HOC1 and I2 are both nonessential. To check these predictions,we suggest shift experiments of the steady-state concentrations of essential species and experiments on the stabilizing or destabilizing effect of constant inflows of essential species. The shift experiments are performed at a stable focus near a Hopf bifurcation. If two Hopf bifurcations are present, where the concentrations of the species differ significantly (as is the case in cross-shaped phase diagrams), the experiments are performed at the stable focus with higher concentration of the species that will turn out to be of type X. A small, constant perturbation (inflow) of each species is added and the effect on the steady-state concentration of the same species is monitored. If the steady-state concentration of the species is higher (lower) with the constant inflow than without, the respective species is said to be + (-) regulated. According to the Citri-Epstein mechanism, which predicts the unstable EC of type 1 to be dominant, Cl02- is predicted to be (-), I- is predicted to be (+), and HOCl, HIO2, and HOI are predicted to be (+). The effect of a constant addition of a species on the stability of the system is first studied at a stable focus close to a Hopf bifurcation. If the system oscillates with the additional inflow of a species, the species is destabilizing (d). If the system remains stable, the species is stabilizing (s). To ensure that the amount of addition is large enough to cross the bifurcation, the procedure is performed on both sides of the bifurcation with increasing amounts of the species added until a change (d or s) in the stability occurs. On starting on theoscillatory sideof the Hopf bifurcation the constant inflow of a species either makes theoscillations cease (s) or leaves the oscillations unchanged (d). The effect of a species must be the same on both sides of the bifurcation. Themechanism predicts ClOz- to be s, I- to be d, and HOCl, HI02, and HOI to be s. Table IV summarizes the predicted behavior of the species by showing the species, its predicted type, and the predicted experimental behavior. Even if not all species can be measured, one may distinguish which type of EC best models the role of species. Furthermore, by means of the experimental procedure given above it is possible to check whether the predicted dominant extreme current shown in Figure l a models the role of the species of the real system satisfactorily. 3. Phase Response Behavior. When perturbations are applied to limit cycles of oscillators, the response of the oscillator is a characteristic phase shift. Phase response curves (PRC) of chemical oscillators and the related phase transitionscurves (PT'C) have become a common tool to investigate the dynamics of

Strasser et al.

2860 The Journal of Physical Chemistry, Vol. 97, No. 12, 1993 chemical reaction mechanisms. Many investigations3k36have dealt with the relaxation oscillationsof the Belousov-Zhabotinski reaction. We apply single-pulse perturbations in calculations to the limit cycles of the Citri-Epstein mechanism9 and suggest from this analysis some further experiments. Phase response studiesprovide more species-specific predictions than do bifurcation studies. It is desirable to perturb as many species as possible, in particular intermediates. PTC are obtained for the mechanism in Table I by calculation of phase shifts caused by single-pulse perturbations applied at different phases of the limit cycle. The fast return to the limit cycle makes such a analysis possible. We choose the rate coefficients as given in Table 11; [C1O2-]0and [I-lOequal 3.187 X l t 3 and 8.563 X 10-3 M, respectively, and we pick ko to be 1.512 X 10-3s-1, For this value of ko the mechanism shows fully developed stable oscillations (see Figure 4). The limit cycle is perturbed by increasing the numerical value of the intermediate HOCl and simultaneously monitoring the response in the Ivariable. We also calculate the phase responses of perturbations with ClOz-. As a reference event for determining phases we choose the very sharp drop of the concentration value of I-. In the following we briefly describe the definitionsof PRC and PTC, and then we focus on the results and interpretations of the calculated phase responses, together with suggestions for possible experiments to check the model predictions. ( a ) Definitions and Terminologies. Let t , denote the time of the nth reference event. Then the period TOof the unperturbed limit cycle is given by

To = t, - t,,

(7)

The reference period used in the calculation is chosen to be the period of the preceding cycle. Assume the time of stimulation, ts,im,is the time at which the pulse perturbation is applied, between t, and rn+l. Then the phase of stimulation is The time of occurrence of the following reference event is generally no longer z , + ~but tm. For t, > z,,+~we obtain a phase delay by perturbing the system, for tm < tn+l we refer to the observation as a phase aduance. The dimensionless entity Acp used as the ordinate of the PRC is defined as

(9) Note that there is an additional minus sign on the right-hand side, unlike the corresponding definition of Acp in ref 36. A PTC is obtained from a PRC by adding cpstim to thevalue of theordinate. Therefore the phase shift plotted on the ordinate of the PTC, defined by d , is

(b) Numerical Results and Suggestions for Experiments. Figure 8, parts a and b, shows the PRC and PTC for HOCl M, respectively. perturbation strengths of 10-4 and 5.0 X For small (large) valuesof cpstim the perturbations delay (advance) the phase of the unperturbed cycle. The transition from phase delay to phase advance occurs according to an abrupt "jump" of the phase shift. Figure 8,parts a and b, show the transitions at qstim = 0.35 and at qstim = 0.25, respectively. Furthermore Figure 8 shows that the absolute values of the phase shifts increase with increasing perturbation strength.36 In Figure 9 the dependence of the phase at which the transition from delay to advance occurs isplottedversus thelogarithmof the perturbation strength. Unlike the results in ref 36 for a mechanism of the BZ reaction, the phase of the transition for the Citri-Epstein mechanism decreases for smaller perturbations. Phase portrait calculations, which are not given here, show that the transitions to phase advance occur during the stage in which the concentration of I- is increasing.

0.6

-0.4

__

1

t

...

0.0

0.2

0.4

0.6

1.0

0.8

pstim

0.4

I

1 *

-0.4t

-064 00

:

:

.

02

:

.

04

: . : 06 08

1

:

10

'Pstim

Figure 8. PRC (curve of triangles) and PTC (Curve of circles) of the Citri-Eptein model for different perturbation strengths of HOCI. The phase shifts Acp and d are plotted versus the phase of stimulation for the PRC and the PTC, respectively. In (a) the perturbation strength equals 1.0 X lV, in (b) it is 5.0 X ko = 1.512 X lo-? [CIO*-]o and [I-]oequal to the values given in Figure 4. Negative Acpcorraponds to a phase delay, positive Acp corresponds to a phase advance. 0

5

1

A A 0.3

A

a

o'2 0.1- 5 . 5

-5.0

-4.5

-4.0

-

1.5

log [perturbation strength]

Figure 9. Dependence of the transition phase from phase delay to phase advance on the logarithm of the perturbation strength of HOCI. The phase of transition decreaseswith decreasing perturbationstrength,unlike the results of ref 36 for models of the BZ reaction.

The interpretation of the phasedelay (phase advance) for small (high) values of cpstim is straightforward: The HOCl perturbations decrease the I- concentration via reaction 5. For small values of cp the perturbations oppose the unperturbed behavior of the limit cycle and cause a delay, whereas for higher values of qstimthe perturbations support the unperturbed behavior and advance the limit cycle. The following considerations explain the unusual dependence of the transition phases from delay to advance on the perturbation strength (see Figure 9): Phase advance is caused by a premature induction of the autocatalysis of HOI and the simultaneous consumption of I-. Using the results of the SNA analysis of section 111, we consider the induction of the autocatalysis as a sudden switch in the dominance of the dynamical elements DE-2 and DE-1: DE-2, consisting of the two inflow reactions of C102- and I-, is dominant during the stage in which I- increases, whereas DE-1 is dominant during the autocatalytic increase of HOI. According to section I11 a switch from an increase of I- to an autocatalytic decrease of I- is successful if

The Journal of Physical Chemistry, Vol. 97, No. 12, I993 2861

The ChloritbIodide Reaction. 1

v1

0.6

t

1

I

t

I

A

0.4

EC 1

0.2

Av 0.0 -0.2 -0.41

A'

A

-0.6 0.0

0.2

0.4

0.6

0.8

1.0

vstim

Figure10. PRC (triangles)and PTC (circles) for perturbationsof [C102-]. The perturbation strength is 1 .O X 10-4. Axes and parameters as given in Figure 8.

the concentrations of HOI and C102- are high enough to turn on the autocatalysis. Simultaneously the concentration of I- may not decrease below a certain value, since I- is crucial to maintaining the autocatalysis via reaction [6]. The perturbations in HOCl provide the necessary increase of the concentration of HOI via reaction [51. Strong perturbations, however, simultaneously decrease [I-] strongly via reaction [5] so that the autocatalysis of the dynamic element DE-1 cannot be turned on. Therefore strong perturbations are only successful at higher values of pstim, where [I-] is higher, whereas weak perturbations turn on the autocatalysis of DE-1 at smaller values of qstim, since [I-] is not decreased as much. Figure 10 shows the calculated PRC and PTC of single pulse perturbations of ClO2-. As before, the response is measured by monitoring the I-concentration. The perturbations, with strength 10-4 M in Figure 10, provide a qualitatively similar PRC as the corresponding HOCl perturbations. The calculations of PRC for various perturbation strengths of C102- predict the same qualitative dependence of the phase of transition from delay to advance on the logarithm of the perturbation strength as has been predicted for the HOCl perturbations (no plot shown): the smaller the perturbation strength, the smaller the value of %tim of the transition. Figure 10 shows, however, that for the same perturbation strength of C102- the transition occurs at a smaller value of vstim ( e t i m = 0.26 in Figure 10, vstim = 0.35 in Figure 8a). To check these predictions, phase response experiments should be performed on the oscillatory side of a Hopf bifurcation, where the system exhibits stable, sustained oscillations monitored by measuring the instantaneous concentration of I-. To define a phase, the fast drop of [I-] is taken as a reference event. The period of the unperturbed limit cycle is determined according to eq 7. Apply HOCl and C102- pulse perturbations by adding alkaline solutions of NaOCl and of ClO2-, respectively, to the stronglyacidicreactionmixture, startingat cpstim = 0 and increasing the value of the phase of stimulation %time The time interval between the previous and the consecutive reference events is monitored, and the phase shift Acp is determined by means of eq 9. The model predicts Acp < 0 for 0 < %tim < & ,:' with cp:?; denoting the phase of the transition from phase delay to phase advance, and Acp > 0 is predicted for qstim > cp:?:. The transition is very abrupt. Near cp::: the value of Acp is predicted to have a minimum (maximum) on the side of the phase delays (phase advances). The value of cp;;: is smaller for pulses of C102-than for pulses of HOCl, with both pulses having the same strength. For smaller perturbationsthe model calculationsshow a smaller absolute value of the phase shifts Acp. Furthermore is predicted to decrease for smaller perturbation strengths of both HOCl and C102-(see Figure 9). By means of Acp and eq 10, the PRC and PTC curves can be determined.

cps:

EC2V

Figure 11. Three-dimensional reaction velocity space spanned by the reaction rate axes ul, u2, and u3 with a two-dimensional steady-state current cone eo.The shaded, closed portion of the originally open cone e, is obtained by intersecting the open cone with the hyperplane Cp, = 1. The one-dimensional intersectionsof the closed portion of @, with the coordinate hyperplanes constitute the two scaled extreme currents ECl and EC2.

A c k n o w l e d ~We thank Dr. Igor Schreiber for illuminating discussions and for making the SNA software package available to us with advice and help. This work was supported in part by the National Science Foundation and the Air Force Office of Scientific Research. Appendix: Essential Ideas of Stoichiometric Network Analysis37 Stoichiometric network analysis (SNA) is an approach to the qualitative analysis of the nonlinear dynamics of networks. SNA was presented by Bruce Clarke in 1980'3 and was used by him for the analysis of various chemical mechanisms.14J2 A chemical mechanism constitutes a set of processes (chemical reactions) with rate laws that remove or produce certain entities (chemical species) in fixed ratios which are given by the stoichiometries of the elementary reactions of the model mechanism. Hence such a mechanism constitutes a stoichiometric network and can be analyzed using SNA. A major objective of SNA is to determine what kind of dynamical behavior is possible for a chosen network. SNA predicts the potential stability or instability of networks without using any rate constants. The range of application of SNA, however, reaches beyond pure stability studies. In the following we briefly describe the essential definitions and ideas of SNA13t36that are used in this study. We assume a network containing n species, Xi,and r reactions, Rj. Xi is the symbolic notation for the species and simultaneously for their concentrations. All reversible reactions are written as two irreversible reactions (see Table I). The net production of species Xi in reaction Rj is denoted by vu, Le., the stoichiometry. All the vu form an n X r matrix Y , called the stoichiometric matrix. Similarly the reaction order of a speciesXi in reaction Rj is given by the element xu of a n X r matrix K, called the kinetic matrix. We further define the reaction ratevector v(X,k) to be a column vector containing the reaction rates of the r reactions:

where kj is the reaction rate constant of thejth reaction and the jth entry of the rate constant vector k. Since kj and Xj are always positive, thevector v(X,k) represents a point in the positiveorthant of the reaction rate space (also called reaction velocity space). Thedifferential equationsthat describethe evolution of a chemical

2862 The Journal of Physical Chemistry, Vol. 97, No. 12, I993 system are

I

dXi/dt = c v i , u j

i = 1,

...,n

(A31

Strasser et al. network contains an unstable feedback cycle and if so where the unstable feedback cycle is located within the network. If the dimension of e, is less than the number of EC of the reaction network, the vector j is not uniquely determined. The procedure given in section 111.2 describes how j can be made unique and how the j i can be interpreted.

j= I

dXldt = vv (A41 Usually at this point the dynamic variables are scaled with respect to their steady state values Xiowhich, however, is not of importance for our considerations. In particular we are interested in the analysis of a chemical system at steady state. Therefore we write 0 = vvo Any solution vector vo is referred to as a "current" of the network. All solution vectors vo form an open subspace @, called 'current cone" or "steady-state cone" which is of dimension ( r - n), if there are no conservation constraints of atoms or molecules. Figure 1 1 shows @, for a system of three reactions and one species. Since all uj are always positive, the current cone is convex. To each current there corresponds a (sub)network with reaction rates UI, ..., u,. For the Jacobian matrix J of a system at steady state we obtain

J = v(diag vo) K'(diag h)

(A61

with hi = l/Xio;t = transpose We use the following method to find the complete set of fundamental currents, called extreme currents: Since@, is convex, any current vector vo can be expressed as a linear combination of vo lying on the boundary of @". Generally the boundaries of e,are facets which are also convex cones but of a lower dimension. The vo on the boundary are simultaneously part of one or more coordinate hyperplanes. Therefore one or more of their uj are zero, which means the vo on the boundary constitute subnetworks that are simpler than the original network since one or more reactions are missing. The given procedure of expressing currents as a linear combination of currents of the boundary can now be repeated with the lower dimensional cones. This leads eventually to the simplest cones, to the one-dimensional edges of e,, called extreme currents. In Figure 1 1 the extreme currents (EC1and EC2)are simply the intersectionsof the two-dimensional @,with the coordinate hyperplanes. From these considerations it follows that any vo E e, can be decomposed into a linear combination of extreme currents with nonnegative coefficients that can be viewed as weighting factors. We can write

vo = Ej where E and j denote the matrix of extreme currents and the nonnegative weightings, respectively. The columns of E are the extreme currents and, as we have seen, the edges of 8,.E is unique up to an arbitrary scaling factor for each extreme current. To each column of E there corresponds an extreme network. The algorithm to calculate the matrix E indicated above is given in more detail in ref 13. Furthermore it is possible to determine the stability of an extreme network with given values of E and j by investigating subdeterminants of the Jacobian of the linearized system.13 It can be determined whether or not a given extreme

References and Notes (1) Kern. D.: Kim. C. H. J . Am. Chem. Soc. 1965.87. 5309. (2j Orbin, Dako, C. E.; DeKepper, P.; Epstein, I. R. J . Am. Chem. SOC.1982, 104, 591 1. (3) Weitz, D. M.; Epstein, I. R. J. Phys. Chem. 1984, 88, 5300. (4) Menzinger, M.; Boukalouch, M.; DeKepper,P.; Boissonade, J.; Roux, J. C.: Saadaoui. H. J. J. Phvs. Chem. 1986. 90. 313. ( 5 ) Luo, Y.;Epstein, I.-R. J . Chem. Phys. 1986, 85. 5733. (6) Ouyang. Q.; Castets, V.; Boissonade, J.; Roux, J. C.; DeKepper, P.; Swinney, H. L. J . Chem. Phys. 1991, 95, 351. (7) Laplante, J. P.; Borkmans, P.; Dewel, G.; Gimenez, M.; Michean. C. J . Phys. Chem. 1987, 91, 3401. (8) Epstein, I. R.; Kustin, K. J. Phys. Chem. 1985.89, 2275. (9) Citri, 0.;Epstein, I. R. J . Phys. Chem. 1987, 91, 6034. (IO) De Kepper, P.; Boissonade, J.; Epstein, 1. R. J . Phys. Chem. 1990, 94, 6525. (11) Olscn, R. J.; Epstein, I. R. J. Chem. Phys. 1991, 94, 3083. (12) Song, W. M.; Kustin, K.; Epstein, I. R. J . Phys. Chem. 1989, 93, 4698. (13) Clarke, B. L. Ado. Chem. Phys. 1980,43, 1-215. (14) Aguda, D. B.; Clarke, B. L. J . Chem. Phys. 1988,89,7428. (15) Aguda, B.; Larter, R.; Clarke, B. J . Chem. Phys. 1989. 90,4168. (16) Stemwedel, J. D.; Ross, J. J . Phys. Chem., in press. (17) The package of programs for solving problems of SNA is kindly

hi.;

provided by Dr. Igor Schreiber. The package comprises the routines SNACONE, SNABETA, SNAEVLS. and SNADCMP. (18) Dcedel, E. J.; Kernevez, J. P.; AUTO Software for Continuation and Bifurcation Problems in Ordinary Differential Equations; AUTO 86 User Manual, Technical Report, Applied Mathematics, California Institute of Technology, 1986. (19) The total flow rate ko used in this study is normalized with respect to the volume of a vessel of 50 cm3. Therefore ko has the units of a reciprocal time. The value of ko constitutes the reciprocal value of the residence time. (20) Lsode: Livermoresolver for ordinary differentialquations, IVPsolver for stiff and nonstiff problems, based on the Gear packages, version of Mar 30, 1987, written by A. C. Hindmarch, Livermore National Laboratory, Livermore, CA. (21) Eiswirth, M.; Freund, A.; Ross, J. Adu. Chem. Phys. 1991,80, 127. (22) Aguda, D. B.; Clarke, B. L. J . Chem. Phys. 1987.87, 3461. (23) Smoes, M. L. J. Theor. Biol. 1976,58, 1. (24) Smoes, M. L. J. Chem. Phys. 1979, 71,4669. (25) Maselko. J. Chem. Phys. 1983, 78, 381-389. (26) Maselko. J. Chem. Phys. 1982, 67, 17-26. (27) Noszticzius, Z.; McCormick, W. D.; Swinney, H. L. J . Phys. Chem. 1989, 93, 2796. (28) Boissonade, J.; DeKepper, P. J . Phys. Chem. 1980,84, 501. (29) Dateo, C.; Orbh, M.; DeKepper, P.; Epstein, I. R. J . Am. Chem. Soc. 1982. 104, 504. (30) Ringland, J. J . Chem. Phys. 1991, 95, 555. (31) See: AUTO 86 User Manual, Technical Report, Applied Mathematics, California Institute of Technology, 1986. (32) See ref 10, section 11, 'Flow (CSTR) Dynamics". (33) Schreiber, 1.; Ross, J., manuscript in preparation. (34) Ruoff, P. J . Phys. Chem. 1984,88, 2851. (35) Ruoff, P.; Noyes, R. M. J . Chem. Phys. 1988,89,6247. (36) Ruoff, P.; Foersterling, H. D.; Gyorgyi, L.; Noyes, R. M. J. Phys. Chem. 1991, 95, 9314. (37) Clarke, B. L. Cell Biophys. 1988, 12, 237. (38) Andronov, A. A.; Vitt, A. A.; Khaikin, S.E. Theory of Oscillalors; Pergamon: Oxford, 1966. (39) Noszticzius, Z.; Stirling, P.; Wittman, M. J . Phys. Chem. 1985.89, 4914. (40) Sheintuch, M.; Luss, D. Chem. Eng. Sci. 1987, 42, 41. (41) Sheintuch, M.; Luss, D. Chem. Eng. Sci. 1987, 42, 233. (42) Noszticzius, Z.; Stirling, P.; Wittmann, M. J . Phys. Chem. 1985,89, 4914. (43) Bar-Eli, K.; Noyes, R. M. J . Chem. Phys. 1987,86, 1927. (44) Brons, M.; Bar-Eli, K. J. Phys. Chem. 1991, 95, 8706. (45) Gaspar, V.; Showalter, K. J. Chem. Phys. 1988,88, 778.