J. Phys. Chem. 1993,97, 6116-6181
6776
FEATURE ARTICLE Toward a Systematic Determination of Complex Reaction Mechanisms Tim Chevalier, Igor Schreiber, and John Ross’ Department of Chemistry, Stanford University. Stanford, California 94305 Received: February 23, I993
Given a complex chemical reaction with an unknown or a partially known mechanism, we examine the possibility of determining the essential parts of the mechanism by a number of experimental methods, all of which rely upon the experimental evaluation of the stationary-state Jacobian matrix elements (JMEs). If the response by all species to a pulse perturbation of one of the species can be measured, then it is possible to determine all of the JMEs. A different experimental method, a concentration shift experiment, relies on measuring the change of steady-state concentrations in a chemical reactor after the inflow of each species has been altered. Another technique employs the induction (or cessation) of oscillations caused by a delayed feedback imposed on an inflow species. If we assume power law kinetics, the connectivity of the reaction network can be determined on the basis of the signs of the JMEs. Once the connectivity of the reaction network has been established, then, with some knowledge of basic chemistry, the essential parts of the mechanism can be constructed. Furthermore, the rate coefficients of the suggested reaction rates can be calculated from the measured JMEs. The method of constructing a mechanism is illustrated with JMEs calculated from a model of the horseradish peroxidase reaction. If the inverse of the Jacobian (the concentration shift matrix) or the phase shifts between pairs of oscillating species can be measured, the species in the oscillatory reaction can be classified and the reaction mechanism can be categorized. If, due to experimental constraints, the entire Jacobian cannot be evaluated, these methods offer a new series of tests that must be passed by any proposed mechanism. Qualitative determination, such as the signs of JMEs, can also be useful.
I. Introduction Oscillatory reactions, multiple stationary states, quasiperiodicity, or chaos occur in complex reaction mechanisms.’-3 Such mechanisms may feature some combination of positive and negative feedback, autocatalysis, and, for ideal reactants (ideal gases or solutions), far-from-equilibrium constraint^.^,^ The formulationof a mechanism for some kinetic systems is a difficult enterprise. This is especially true for reaction sequences involving numerous species and many elementary steps, as is frequently found for oscillatory reactions. In this paper, we focus on determining essential parts of a mechanism for an oscillating homogeneous isothermal reaction system. We also show how the reaction speciesin the reaction can be classified and the mechanism categorized using information derived from the inverse of the Jacobian matrix or from the phase shifts for the oscillating species in the reaction. This paper is a further contributiontoa systematic approach leading toward the determination of complex reaction mechanism^.^ Suppose we have a mechanism involving n speciesfor a reaction occurring in a reaction cell. The evolution equations for all of the species are dXi/dt = Fl(XI,..,Xn), i = 1, ..., n
(1)
where Fiis, in general, a nonlinear function of the concentrations X1,...,X,. For a stationary state A? = 4, ...,xl, we have
Fl(x”) = 0, i = 1, ..., n
(2)
An experimentallyimportant way of keeping the system far from equilibrium is to use a CSTR (continuous-flow stirred tank
0022-3654/93/2097-6776$04.00/0
reactor).*O In this case, Fl can be expressed as
,
where Ri represents the change of species i due to chemical reactions, ko is the reciprocal residence time, and is the feed concentration for the ith species. We may expand eq 1 in a Taylor series about a reference solution x7 = q ( t ) , ...,x(t), obtaining to first order
Here, SXf = Xl- X: and the partial derivative Jij = aFi/aXj is a Jacobianmatrix element (JME). A particular choiceof interest f o r m is a stationary s t a t e x . The sign of Jij determineswhether theconcentrationofspeciesi willbeincreased(JU> 0) or decreased ( J f j< 0) by adding a small amount of Xi. The numerical values of the JMEs are related to the mechanism of the reaction through the stoichiometric and rate coefficients and the concentrations of the species. Thus, measurements of the JMEs provide information on the mechanism of the reaction being investigated. Conversely, the calculated predictions of the JMEs from a stipulated reaction mechanism provide a substantial set of predictions (more extensive than a concentration time series) which need to be confirmed by experiments prior to acceptance of that mechanism. Tyson6 has discussed relations among the JMEs that need to be satisfied so that either positive or negative feedback loops, or both, produce an instability. A classification schemefor oscillating chemical systems, based upon theinterplay of positive and negative feedback loops, has been suggested by Luo and Epstein;2another 0 1993 American Chemical Society
Feature Article one, based on stoichiometricnetwork analysis? was proposed by Eiswirth, Freund, and RossS3 We present in this article several methods for the partial determination and classification of oscillatory reaction mechanisms. In section 11, we describe several new methods for determining the JMEs from measurements and cite briefly, for completeness, prior methods reported in the literature. In most of the methods all of the species in the reactive system must be measurable if all of the JMEs are to be determined. This is hardly ever attainable; however, as we shall show, subsets of the Jacobian matrix can provide important information. Qualitative determinations (the signs of the JMEs) can also be very useful. One method, which involvesthe use of a delayed feedback, requires the measurement of one species only. In section I11 we choose a reaction mechanism and evaluate the Jacobian matrix near a Hopf bifurcation. We then use the same reaction mechanism as a “black box” experiment to evaluate the JMEs by the various methods presented,consider the likely uncertainties encountered in each method, and compare them to the analytically obtained Jacobian. Then, in section IV, we show how this information may be used to construct a reaction mechanism, to determine reaction rate coefficients, to identify the roles of the chemical species in the reaction, and, finally, to categorizethe mechanism.3 The determination of the JMEs leads to a skeletal reaction mechanism with the roles of key chemical species in that mechanism identified; furthermore, it provides tests to be passed by a proposed mechanism. The classification of the mechanism may be based upon the elements of the inverse of the Jacobian matrix or on simple measurements of the mutual phase shifts of the oscillating species.
11. Methods for Determining the JMES In the following, we present three methods for determining the JMEs: pulse perturbations (section ILA), concentration shift experiments (section ILB), and a method based upon delayed feedback experiments (section 1I.C). In each of these cases, we assume that a chemical reaction, with a complex reaction mechanism, is kept far from equilibrium in either a stationary state (sections II.A.l, II.A.2, II.B, and 1I.C) or on a limit cycle near a Hopf bifurcation (sections II.A.3 and 1I.C). The system can be maintained in such a state by running the reaction in a CSTR. The flow-through arrangement is not essential; some species, such as enzymes, for example, may be confined to the reactor vessel. The flow terms are included in the Fls, or noted specifically [see eq 31. A. Methods Using a Pulse Perturbation. There are at least three different types of pulse perturbation experiments. The first two methods discussed (sections II.A.1 and 2) assume that the kinetic system is in a stationary state which is perturbed upon the addition of species kat a given instant to. In this situation, provided that the perturbation is small so that the first-order approximation by linearized kinetic equations is sufficient, eq 4 describes the relaxation back to thestationary state. The third method (section II.A.3) may be applied to limit cycle oscillations located near a Hopf bifurcation. In this case, a pulse can temporarily quench the oscillations by bringing the system close to the unstable stationary state; this serves as the basis for employing the linear equations (4). In order to solve eq 4 for the JMEs, we need to measure the species that significantly affect the dynamics of the system; products that do not reenter any chemical reaction and species that do not show any temporal variation, because of an imposed constraint of constant concentration (a pool condition), need not be monitored. The rest of the species, which we call essential species,3 must be accessible to measurement, if we want to determine a complete Jacobian. However, in practical applications, the investigated mechanism is almost never a black box, so that knowledge of a part of the mechanism may be sufficient to avoid measurements of all essential species. In any case, it is best to measure as many species as possible.
The Journal of Physical Chemistry, Vol. 97, No. 26, 1993 6777 1 . Method of Initial Slopes. Bar-Eli and Geiselera have considered a perturbation of a kinetic system in a stationary state of sufficiently small amplitude such that eq 4 is valid; only a single species, say the kth one, is perturbed by injection into the CSTR. The decay of that perturbation is assumed to be given by
Comparing with eq 4, we can see that eq 5 is exactly true only at the initial time of addition of species k,since only then are all other SXJ’szero. The initial slope of the temporal evolution of the species directly gives the kth column of the Jacobian. In order to obtain the full Jacobian, all of the species need to be perturbed and measured. If np species are used as perturbants and n, species are measurable, then we obtain a npX n,submatrix of J. A reliable extrapolation of the initial slope may require very fast measurement techniques. For the system studied by Bar-Eli and Geiseler-the oxidation, in the presence of B r , of cerous or manganous ions by bromate in H2S04-this method requires sampling periods on the scale of micro- to milliseconds to be useful. 2. Method of Multilinear Least-Squares Fit. Let us assume that we can measure the temporal evolution of the concentrations Xl, ...,Xnof all dynamical species in the form of a time series. We begin at the steady state (4,...&). After applying a perturbation, by injecting a small amount of a selected species into the reaction vessel, we obtain a vector time series ,...tx’,),1 = 1, ...,N . This time series can be conveniently put into a matrix U such that Ut1 = We assume that the number of measured data points N is large, N >> n, and that the sampling rate is dense enough to cover the fastest time scale of the system. Without assuming any particular reaction mechanism, we can solvefor the Jacobian matrix elements in (4) by applying regression analysis. Instead of theexact solution (SXl(t),...,SXn(t)),weuse the matrix 6U obtained from U by subtracting the steady state values, Sui,= Ut1 - q. The time derivatives, dSXi(t)/dt, are approximated by a matrix dbU/dt which is obtained from 6U;the easiest, but not necessarily the best, way to obtain the derivatives is by numerical differencing. After the substitution of 6U for SX(t), eq 4 reads
(4
4.
d6U/dr = J6U where the unknown is the Jacobian matrix J. As the number of measured data points (columns of 6U) is large, this linear system is overdetermined. Due to the finite precision of the measurements, we cannot expect a (classical) solution to (6) to exist. However, a least-squaressolution that minimizes the norm of the difference d6U/dt - J6U is defined uniquely. The problem, therefore, is reduced to a standard multilinear least-squares fit problem. The system of equations may be solved by singular value decomposition.9 This procedure solves for one row of J at a time with the use of a pseudoinverse of UT. The accuracy of the estimates of the Jacobian matrix elements obtained with the outlined method is dependent both on the precision of the measurements and on the sampling rate. With a limited precision of measurements, the numerical differencing procedure may generate large errors that would invalidate the solution of eq 6. This problem can be overcome by using an integral representation of eq 4, with the solution expressed in terms of the eigenvectors u1, ..., un and eigenvalues AI, ..., A,, of the Jacobian matrix J n
(7) Here, the unknown parameters WJ = a p i and A, must be obtained by nonlinear regression. A numerical procedure based on the least-squares estimation, such as the Levenberg-Marquardt
6778 The Journal of Physical Chemistry, Vol. 97, No. 26, 199'3 method? can be used for this purpose. Drawbacks of this approach are that the nonlinear regression does not necessarily provide a unique solution and that finding a good fit can be quite laborious. After fitting the measured data N)to (7), the time derivative can be determined through the analytical differentiation of eq 7. To recover the Jacobian, the estimated parameters in (7) can be used in two ways. One method of recovering the Jacobian depends upon calculatingd6Uldt (via the analyticaldifferentiation of (7) and the definition of 6U), followed by the use of eq 6 . Or, the Jacobian can be directly evaluated as
J = WAW'
(8) where W = ( w A and A is a diagonal matrix containing the eigenvalues. The sampling rate is also an important factor, particularly in systems possessing either very fast time scales or vastly different time scales. The problem of vastly different time scales is a major obstruction in using the outlined method. In such cases, after the perturbation, the system quickly relaxes to a submanifold (dimension less than n) of the concentration space. Within the neighborhood of the steady state governed by the linearized equations, this submanifold is a linear subspace spanned by a subset of eigenvectors of J. If the relaxation of the transient to this linear subspace is too fast to be measured accurately, the matrix 6U is almost singular and eq 6 cannot be applied; see section I11 for further discussion of this point. For special choices of initial perturbations, only some of the JMEs may be obtainable, since not all of the species need to be affected by a given perturbation. The pulse perturbation method requires the experimentalist to measure all specieswhile perturbing the reaction with one species. 3. Quench Experiments. Each of the methods discussed so far requires perturbing a stable stationary state. It is also possible to determine the JMEs by perturbing a stable limit cycle. We assume that the limit cycle oscillations appear via a Hopf bifurcation, implying that there is an unstable stationary state nearby. Close to the Hopf bifurcation, there exist both a twodimensional unstable manifold and a dimension-two stable manifold of A?. If we manage to perturb the limit cycle so that the shifted system hits the stable manifold of 2, then the oscillations may be quenched.*O Due to the ever-present fluctuations in the system, the experimental quenching is always temporary; eventually, the limit cycle oscillations are restored. Since eq 4 is valid only in a small neighborhood of A?, we need to take that portion of measured data which falls within a small neighborhood of A?. Once we have selected the appropriate portion of the measured data, the calculation of J proceeds in the same fashion as described in section II.A.2. In any case, an almost exact quenching is required to find the stationary values 4, ..., Y,, needed for the evaluation of 6U. The quenching of oscillations has been used10 with the idea of getting as much information from an experiment as possible, provided that all species can be added but only one species can be measured. The quenching is described in terms of both the amount of species added and the phase of the oscillation during the addition of perturbant. If the limit cycle is very close to the Hopf bifurcation, the oscillations are nearly planar. One may attempt to use eq 4 to calculate the pair of complex conjugate right eigenvectors associated with the pair of pure imaginary eigenvalues of J. These eigenvectors span the plane containing the limit cycle. With quenching data for all dynamical species, it is possible to calculate a pair of left eigenvectors of J. It is necessary to measure n - 2 oscillating species (rather than 1) if a pair of the right eigenvectors of J is desired. A complete determination of the Jacobian matrix requires knowing all of the eigenvalues and eigenvectors. B. ConcentrationShift Measurements. Useful measurements for inquiries into the categorization of oscillatory reactions of a particular system are the so-called shift experiments.3 Consider
Chevalier et al. the CSTR arrangement described by eqs 1 and 3. Suppose that at a stationary state A? we alter the inflow concentration of species j by 6x0; after a transient, a new stable stationary state is assumed to set in. The change in the stationary concentration of species i in the CSTR relative to the perturbation, 6x/6$, can be well approximated by *Id$, obtained by differentiating eq 2 with respect to 3,
From this follows (in matrix notation)
A X = W / d p = -k,(dF/dX)-'
(10)
If the input concentration shifts a$ are sufficiently small, then the shift matrix AX is proportional to the inverse of the Jacobian matrix J. The JMEs can be obtained consequently by inverting the shift matrix. Moreover, AX by itself contains useful informationconcening the categorizationof oscillatory reactions; see section IV. The shift experiments can also be performed in non-CSTRsystems,provided that thereare externallycontrollable zero-order terms in eq 1. For example, the shift on input can be realized by turning on a slow infusion of a species. The relation of AX to J is based on the assumption that the Jacobian is not singular. Except for marginal cases, this is always fulfilled for the CSTR but may be violated for other systems. In particular, singularity is present in systems with conservation constraints. Here, the shift perturbations can be applied only to the species which are not bound by the constraint(s) or else there is no stationary state. The Jacobian of eq 1 can be made regular by the standard procedure of eliminating some species with the use of the constraints. In these situations, the shift matrix is proportional to the inverse Jacobian for the reduced system. C. Delay Methods. A delayed feedback in an experimental chemical system can lead to new dynamical behavior which is not present without the feedback.11-13 For example, using a delayed feedback, it is possible to introduce oscillations into a system that doesnot exhibit oscillationswithout the feedback. Theappearance ofoscillationsin thesystem with feedbackcanbeused todetermine quantitatively the elements of the Jacobian of the steady-state system. I , Stability in Delay Systems. As is the case with ordinary differential equations, the stability properties of solutions to a differential delay equation can be determined by linearizing the delay equation about the steady states of the system and then examining the solutions of the characteristic equation for the linearized system. For a number of differential delay systems (including chemical systems with power law kinetics and one imposed delayed feedback), the characteristic equation has the simple form
P(z) + Q(z)e-" = 0 (11) where P and Q are polynomials of degree n and m, respectively, with real coefficients and T is the imposed delay. As the delay T is increased, the stability properties of the delayed system may change. For characteristic equations with n > m (the case for chemical systems with an imposed delayed feedback), an increase in 7 may lead to a series of stability changes or stability switches as the zeroes of eq 11 cross the imaginary axis. Cooke and van den Driessche14have outlined a method for determining the stability properties of solutions to the characteristic eq 11 when n > m. With the provisions (1)-P ( z ) and Q(z) have no common imaginary zero; (2) P(-iy) = P(iy),Q(-iy) = Q(iy) for real y ; (3) zero is not a root of (1 1); (4) eq 11 has at most a finite number of roots in the right half of the complex plane for T = 0; and ( 5 ) FCy) IP(iy)l2- lQ(iy)12 has a finite number of real zeroes for real y, they find that (a) if FCy) = 0 has no positive roots, the stability properties of (1 1) are not
-
The Journal of Physical Chemistry, Vol. 97, No. 26, 1993 6779
Feature Article changed by the delay; otherwise, (b) as the delay increases, a finite number of stability changes may occur. If a real positive root of FCy) is found, then (1 1) holds if and only if
+ (-l)"lko"ltr(Jp) + (-1)"2ko"2tr2(Jp) + ... + k;tr,(Jp) + (-l)kotr,+l(Jp) + tr,,(Jp) (19)
p,, = (-l)"k,"
where P&), PI^), QRC~), and Q1Cy) are, respectively, the real and imaginary part of P(iy) and the real and imaginary part of Q(~Y). 2. Experimental Determination of the Jacobian in Delay Systems. One can use a delayed feedbackto provide quantitative information about a chemical system. Any generaln-dimensional isothermal chemical system reacting in a CSTR can be represented by eq 3. Initially, we assume that only one species in the reaction (for example XI)can be monitored. The delayed feedback is introduced into the system by specifying that at least one of the inlet concentrations $ ( t ) is a function of the concentration of Xl(t-7) at an earlier time t - 7. That is,
$(t) =f(X1(W) (13) With one delay, the characteristic equation for the general iystem represented by eq 3 can be expressed as
=O
:I41
In (19), the trace of the Jacobian J p is represented by tr(Jp). The remaining traces are defined as the sum of the determinants of all (i)matrices of dimension k X k which are formed by the rows and columns of J p such that they intersect at any k-tuple of diagonal elements of Jp. Note that tr,, for an n X n matrix is equal to the determinant. Likewise, the coefficients of Q are also defined in terms of the traces of JQ: 40=
4-1
- (-1)"'k0"'
1
+ (-l)"2ko"2tr(JQ) +
+ (-1)'$)tr&(JQ) + tr,,(JQ) (20) In (20), trk(JQ) is the kth trace of J Q There are a maximum of 2n - 1 unknowns to be determined (the n unknowns tr(Jp) through tr,,(Jp) and the n - 1 unknowns tr(JQ) through tr,l(JQ)). By finding n Hopf bifurcations with n different feedback functions, the unknowns can be determined through the use of eq 12. Once the traces tr(J), ...,tr,,(J), tr(Jp), ..., trWl(JQ)are determined, the eigenvalues for the system can also be determined since
or, equivalently, IJp - I(k0
+ Z)I + ko
e-"lJQ - I(k0
+ Z)I = 0
(15)
where I is the identity matrix, J p is the full Jacobian (without any delay terms), and JQis the reduced matrix of dimension (n - 1) x (n - 1) formed from J p by the removal of the first row and column. Upon being expanded, it is apparent that eq 15 can be expressed in the same form as eq 11. We are interested in the situations in which the chemical system with the delayed feedback exhibits oscillations. If a Hopf bifurcation is to occur, eq 15 must possess an imaginary root. As discussed above, z = iy is a root of eq 15 if and only if y is a positive root of
FoI) = IP(~Y)I~ -lQ
