8556
J. Phys. Chem. 1996, 100, 8556-8566
Categorization of Some Oscillatory Enzymatic Reactions Igor Schreiber Department of Chemical Engineering, Prague Institute of Chemical Technology, Prague, Czech Republic
Yu-Fen Hung and John Ross* Department of Chemistry, Stanford UniVersity, Stanford, California 94305 ReceiVed: September 26, 1995; In Final Form: February 8, 1996X
We investigate the categorization of two or more proposed reaction mechanisms for each of the following oscillatory enzymatic reactions: (1) the peroxidase-oxidase reaction; (2) glycolytic oscillations; (3) oscillations of cyclic AMP in slime mold cells; (4) enzymatic pH oscillations; (5) calcium spiking in cytosol. We use prior work in stoichiometric network analysis and categorization of oscillatory reactions to identify in each proposed reaction mechanism essential and nonessential species, the specific role of each essential species, the connectivity of the essential species, including the identification of the reactions leading to oscillatory instabilities, and the category. For each model, we predict the results of several experiments including relative amplitudes, quench amplitudes, phase shifts, and sign symbolic concentration shifts and compare them with those from available experiments. These and several other experiments such as bifurcation analysis, phase response curves, entrainment experiments, qualitative and quantitative pulsed species response, delay experiments, and external periodic perturbation provide stringent tests of proposed reaction mechanisms, and appropriate ones are suggested to discriminate among competing mechanisms for a given reaction. We find the necessity for introducing a new subcategory in our categorization of oscillatory reactions. The details of the analysis for reactions 2-5 and of the methods determining the categories are given in the Supporting Information.
1. Introduction Chemical reaction systems, which may have oscillatory kinetics of multiple steady states far from equilibrium, usually have complex reaction mechanisms. In prior work1 we have proposed a classification of oscillatory reactions. For reaction mechanisms that contain only one source of instability (based on one source of autocatalysis), we found that we could fit some 25 reaction mechanisms, inorganic and organic, into four defined categories. For each mechanism we first identified essential and nonessential species according to operational definitions. We then used stoichiometric network analysis (SNA)2 to classify reactions depending on their basic unstable subnetworks and dominant feedback. In each of the four categories each essential species is assigned a specific role. Several types of experiments were suggested for the assignment of a given reaction system to one of the four suggested categories, which describe the role and connectivity of the essential species. In further complementary studies3,4 we examined the possibility of determining at least the essential parts of a complex reaction mechanism by various experimental methods, many related to measuring the Jacobian matrix elements of a stationary state near a Hopf bifurcation to oscillations. From such measurements we showed that for power law kinetics the connectivity of the reaction mechanism can be determined, that is, the flow and transformation of essential species in the reaction network. The articles3,4 detail the connection among signs of the Jacobian matrix elements as determined by different types of measurements and the four different elements of categorization. An example (a model of the peroxidase-oxidase reaction4) is used to show that the theoretical analysis provides the basis for determining essential parts of the reaction mechanism of an * To whom correspondence should be addressed. X Abstract published in AdVance ACS Abstracts, April 15, 1996.
S0022-3654(95)02853-X CCC: $12.00
oscillatory reaction from cited types of measurements. Other methods, such as bifurcation analysis,5 measurements of phase shifts between pairs of oscillating species, phase response curves,6 quenching experiments,7 entrainment experiments,8 and determination of dominant subnetworks in the stoichiometric network analysis of a reaction mechanism, are all useful in forming a strategy for the difficult task of establishing a mechanism and verifying a proposed mechanism. A number of the experiments suggested in refs 3 and 4 have been carried out in two systems: the chlorite-iodide reaction9 and the horseradish peroxidase reaction (oxidation of NADH).10 In each, substantial information on the reaction mechanism was deduced from these designed experiments. For the chloriteiodide reaction, the detailed Citri-Epstein mechanism was confirmed with one exception: the identification of the type of two of the nonessential species, I2 and Cl-. For the horseradish peroxidase reaction, the experimental data suggest that the oscillatory system belongs in category 1C, most likely 1CW. This provides a criterion for proposed models of the horseradish peroxidase system. In the present work we focus on the categorization of a few enzymatic oscillatory reactions with some of the methods used in prior work3,4 and suggest how they might be useful for determining the appropriate reaction mechanism for describing the system. Thus, the purpose of this study is a contribution to a systematic approach to the testing and formulation of reaction mechanisms of enzymatic reactions. For other contributions to this systematic approach see ref 11. We study in some detail the following reactions: (1) the peroxidase-oxidase reaction, (2) glycolytic oscillations, (3) oscillations of cyclic AMP in the slime mold Dictyostelium discoideum, (4) enzymatic pH oscillators, (5) and calcium spiking in cytosol. For each of these reactions we analyze proposed mechanisms and use prior work on stoichiometric © 1996 American Chemical Society
Oscillatory Enzymatic Reactions
J. Phys. Chem., Vol. 100, No. 20, 1996 8557 symb TABLE 1: Symbolic Phase Shifts ∆φi,j and Sign Symbolic Concentration Shifts Sign (Dxs,i/Dx0,j) of Prototypes of the Categories 1B, 1CX, 2B, and 2Ca
Category 1B i
j symb ∆φi,j
X Y Z X Y Z
sign(∂xs,i/∂x0,j)
X
Y
Z
I + + + + +
I + -
+ I + -
Category 1CX i
j symb ∆φi,j
X Y Z X Y Z
sign(∂xs,i/∂x0,j)
X
Y
Z
I A + -
A I + + +
+ I + -
Category 2B i X Z X Z
j symb ∆φi,j
sign(∂xs,i/∂x0,j)
X
Z
I + + +
I -
Category 2C i X Z X Z
j symb ∆φi,j
sign(∂xs,i/∂x0,j)
X
Z
I + -
+ I + -
a Symbols I, A, +, and - represent the phase shifts corresond to in-phase, antiphase, advanced, and delayed oscillations, respectively, of species j, with respect to species i. Symbols + and - represent concentration shifts correspond to an increase and decrease, respectively, in steady state concentration of species i upon an increase in inflow concentration of species j. Notice that the shifts for corresponding species in 1B and 2B (1CX and 2C, respectively) are the same. The corresponding species in category 1CW have the same shift behavior as in 1CX, and the additional species W has the same shift behavior as species X.
Figure 1. List of reactions and corresponding network diagrams of prototypes of basic categories and subcategories.
network analysis2 and categorization methods.1,3,4 A brief review of the terminology used throughout this paper and the basis of these techniques are given in the Appendix. The details of the analysis for reactions 2-5 and the Appendix are given in Supporting Information. The analysis leads, in each case of a proposed reaction mechanism, to the identification of essential and nonessential species, the identification of specific roles of essential species, and the identification of the category of oscillator. Once this is done, then the suggested experiments are available to test the predictions of the analysis. These experiments, which include relative phase shifts, concentration shift regulation, concentration shift destabilization, qualitative and quantitative pulsed species response, delayed experiments, quenching, and external periodic perturbation, are more stringent tests of the proposed reaction mechanisms than measurements of the time series of concentration of chemical species.4 Furthermore, from these measurements, crucial features of the reaction mechanism can be deduced.9,10 2. Methods for Determining the Categories In our previous work1,3 we proposed a scheme of classification of known chemical oscillators into two main categories: 1
and 2. Within category 1 oscillators, two distinct subcategories 1B and 1C are found, and the latter is further divided into 1CX and 1CW subcategories. Each of the four basic categories is represented by a prototypical skeleton reaction mechanism. In ref 1 we found no mechanisms that necessitated a division of category 2. In this work, however, we do find mechanisms that require the splitting of category 2 into subcategories 2B and 2C. The earlier category 2 becomes 2C, and a new category 2B is introduced. Examples of mechanisms falling within category 2B are provided by a model of glycolysis and by several theoretical models including a Volterra-Lotka scheme (see section 4). Therefore, in this work we extend the list of categories and give prototypes for each of them in Figure 1. Each prototype is described by a set of chemical reactions in the standard notation of stoichiometry and graphically via network diagrams. The reactions involve four types of essential species: X (autocatalytic), Y (exit), Z (feedback) and W (recovery); see Appendix for further terminology and details. Basic indicators of the role of essential species and of the categories used throughout this paper are the qualitative phase and concentration shifts. For reference the qualitative shift matrices associated with each of the prototypes are shown in Table 1. It should be emphasized that even if prototypes may be chosen in a number of ways (for instance, an exit feedback
8558 J. Phys. Chem., Vol. 100, No. 20, 1996 instead of tangent feedback can exist in 1C), the shift behavior is unchanged. Our approach to the classification of a reaction mechanism is based on the following computer-assisted analysis. (i) Given the stoichiometric matrix ν of a mechanism, determine the extreme subnetworks e1, ‚‚‚, eN of the mechanism and their combinations that form subnetworks still simpler than the entire mechanism. This is done by determining the edges and faces of the stationary state cone that represents the set of all reaction rate vectors at steady state. (ii) Check the stability of these subnetworks in terms of principal minors βk of a stability matrix V associated with the Jacobian matrix J (see eqs A4 and A5 in Appendix). Each negative βk indicates a k-tuple of species that may bring about an instability of a steady state of a subnetwork containing those k species, provided that their steady state values are small enough compared to others. In this way, unstable subnetworks that are the dominant reaction pathways for bifurcations to oscillations or multiple stationary states are indicated. Steps i and ii do not require an explicit knowledge of stationary concentrations (xs1, ‚‚‚, xsn) if power law kinetics is assumed or if the effective reaction orders are known. Partial (and in simple mechanisms even complete) classification can be done already at this point; we find that, given an unstable (sub)network, all species indicated by a negative βk are likely candidates for essential species of type X, Y, or W. Conversely, all type X, Y, or W species are indicated by a negative βk. Thus, the only essential species that need not be present in low concentrations is of type Z. Furthermore, if the autocatalytic cycle can be found by inspection of the network diagram, then type X species can be separated from type Y and W species, and type Y species can be found as the species that enter an exit reaction from the autocatalytic cycle. (iii) Generate a steady state reaction rate vector (current) Vs N R e as a linear combination of all extreme subnetworks ) ∑l)1 l l so that a particular unstable subnetwork eL selected from those found in step ii is dominant (by choosing RL large compared to others). Choose steady state concentrations as suggested by the corresponding βk (values of rate coefficients are then implied) and calculate the Jacobian and its eigenvalues. (iv) Search for Hopf bifurcation conditions, i.e., for a pair of pure imaginary eigenvalues (iω of J, by varying either the Rl’s (but still preserving the dominance of the chosen unstable subnetwork) or the steady state concentrations (xs1, ‚‚‚, xsn). This step can be done either interactively, by varying the conditions based on previous guesses, or automatically by instructing the computer to search systematically through a selected region of parameter space. The first approach gives considerable insight into the relative importance of different pathways in the reaction network. It usually leads to a Hopf bifurcation in a few iterated guesses and has been given preference in this work. (v) At the point of a Hopf bifurcation, the (complex) eigenvectors of J associated with the pure imaginary eigenvalues (iω determine the amplitudes and phases of the oscillating species. The eigenvectors of the transpose of J associated with (iω determine the quench vectors, i.e., a concentration increment (quench amplitude) and a phase for each species such that the addition of a species, increasing its concentration by a corresponding increment delivered at a corresponding phase, quenches the oscillations.7 The inverse of the Jacobian provides the concentration shift matrix ∆x (see eq A7). However, J is singular for most models of enzymatic systems because usually there are conservation constraints. Consequently, the concentration shift perturbations cannot be applied to the species involved in the conservation constraint and the shift matrix is not
Schreiber et al. complete (see eq A8). Likewise, the quench perturbations by addition of constrained species violate the constraint. Thus, the quench amplitudes corresponding to such species should be used with caution when determining the category. (vi) Separate the essential and nonessential species, based on their relative amplitudes of the oscillations and the quench amplitudes calculated from the Jacobian, and determine the role of the essential species and the category of the mechanism from the phase shifts, the concentration shifts, or other methods.1,3 Since the concentration shift matrix ∆x is incomplete for mechanisms with constraints, the phase shift relations among the species are the primary means of the classification. Nonetheless, we found it very helpful to perform additional tests with hypothetical skeleton mechanisms in which one or more species, possibly including those bound by constraint(s), are fixed at their stationary values. Every fixed species reduces the number of dynamical variables (and also the dimension of J) by one. Fixing at least one constrained species per each constraint makes J regular at the cost of violating the constraint(s). As a result, the concentration shift matrix for the reduced model is complete and provides an additional test for classification of essential species and categorization of the full mechanism. Moreover, this is an additional test for determining nonessential species.1 We also find that quench amplitudes calculated for species involved in a constraint help to determine whether the species is essential or not, even if the quenching violates the constraint. If a mechanism has more than one independent unstable subnetwork (associated with edges and faces), the outlined procedure is repeated from step iii with a different unstable subnetwork selected as the dominant one. It is expected that the set of essential species, or their role or even the category, may change when different unstable subnetworks are chosen as the dominant ones. In general, we need to examine all qualitatively different combinations of subnetworks giving rise to oscillations and compare them to available experimental data. Since most of the model mechanisms to be discussed are not detailed and relevant experimental data are scarce, we focus on a qualitative rather than a quantitative comparison with experimental results. For the purpose of the categorization, it is sufficient to find an arbitrary Hopf bifurcation for each qualitatively different situation in a model mechanism. If one of the cases qualitatively matches experimental data, the search for a Hopf bifurcation can be formulated as an optimization procedure and used to match the model with experimentally determined characteristics of the emerging oscillations. For example, one can measure amplitudes and phase shifts of the oscillating species near the Hopf bifurcation, or quench amplitudes and quench phases, and estimate the set of reaction coefficients. So¨rensen and Hynne7 used quench data from the Belousov-Zhabotinskii reaction to estimate reaction coefficients in the Oregonator type of model. In the following sections we apply the outlined analysis to model mechanisms of some oscillatory enzymatic reactions that have been treated both experimentally and theoretically in some detail in the literature. The mechanistic source(s) of oscillatory dynamics in each model is(are) examined and interpreted within the framework of the classification approach. The consistency of the models with available experimental data is discussed, and experiments are suggested for testing several predictions. The section on peroxidase-oxidase reaction is treated in more detail than the rest because more relevant experimental data is available.10 Also, for clarity we include in this section detailed information on the classification procedure, which is made more concise in later sections.
Oscillatory Enzymatic Reactions
J. Phys. Chem., Vol. 100, No. 20, 1996 8559
Figure 2. Network diagram of model A for the horseradish peroxidase (HRP) reaction. The unstable subnetwork here and in all following figures is shown in thick lines.
3. Peroxidase-Oxidase System Experimental observations of the dynamics of the peroxidase-oxidase (PO) system including bistability and periodic, quasiperiodic, and chaotic oscillations.12-21 Several mechanisms have been put forward. (1) The Yokota and Yamazaki model12,13,22 is the first attempt to describe comprehensively the reaction based on the authors’ experimental findings. (2) Another model, called model A by Aguda and Clarke,23-25 is based on a systematic selection of some of the reactions that have been suggested by Yokota and Yamazaki,22 Olsen and Degn,26 and others. (3) A very similar model is due to Fed’kina, Ataullakhanov, and Bronnikova.27 (4) Model C28 is an extension of model A. Recently, several more elaborate extensions of the PO mechanism have been suggested.29-31 In addition, there are several abstract models. The Degn-Olsen-Perram (DOP) model32 and a very similar model by Olsen33 are based on a theoretical consideration about branched chain reactions. The model by Alexandre and Dunford (AD)34 is a highly skeletonized version of the general mechanism considered by Yokota and Yamazaki. 3.1. Model A. The simplest of the realistic mechanisms for the PO system is model A with the following set of reaction steps (see also the network diagram in Figure 2):
per3+ + H2O2 f coI
(PO1)
coI + NADH f coII + NAD•
(PO2)
coII + NADH f per3+ + NAD•
(PO3)
coIII + NAD• f coI + NAD+
(PO4)
per3+ + O2- f coIII
(PO5)
NAD• + O2 f NAD+ + O2-
(PO6)
NADH + O2- + H+ f H2O2 + NAD•
(PO7)
2NAD• f (NAD)2
(PO8)
f O 2, O2 f
(PO9) (PO10)
NADH is assumed to be in surplus and its concentration effectively constant; H+ is buffered; NAD+ and (NAD)2 are products. All these species are nonessential. The system is
externally supplied with oxygen via interface mass transfer from the air (pseudoreactions PO9, PO10). No other species are supplied or removed, and products are assumed not to interfere, even when accumulating. The native form of the enzyme (horseradish peroxidase), per3+, is oxidized by the superoxide ion radical O2- to form oxyperoxidase (coIII), which is successively reduced through enzymatic forms known as coI and coII, respectively, back to the native form. The first step in the reduction (PO4) consumes the radical NAD•, while the second and third steps (PO2, PO3) consume NADH and produce two NAD•. Hence, there is a net stoichiometric production of NAD• in an autocatalytic double-cycle: NAD• f coI f NAD• and NAD• f coI f coII f NAD•. This turns out to be the major cycle participating in oscillations. The cycle is driven directly by NADH via tangent reactions PO2 and PO3, and indirectly by O2, through reactions PO6, PO5, and PO4. Hydrogen peroxide is formed by the entrance reaction PO7, and H2O2 is a competitive substrate to the superoxide ion for the native enzyme by which it is oxidized to coI through reaction PO1. This particular reaction is one of the most thoroughly studied reactions in the PO system (see refs 22 and 26 and references therein). The network has 10 (pseudo) reactions PO1-PO10 and 8 dynamical species: per3+, coI, coII, coIII, H2O2, NAD•, O2-, and O2. Our calculations (see Appendix for methodology) indicate a three-dimensional cone and one conservation constraint, [per3+] + [coI] + [coII] + [coIII] ) const. The three extreme subnetworks involve the following reactions: (i) PO9, PO10, (ii) PO2-PO6, PO9, (iii) PO1-PO3, PO6-PO9. Denote e1, e2, and e3 as the corresponding extreme currents. If we assume mass action kinetics, then we find from the calculated βk’s that while the stationary states of e1 and e3 are always stable, the principal minor β4 associated with the combination of the species coI, coII, NAD•, and O2 and two β5’s involving the previous four species plus per3+ and O2-, of e2 is negative and thus the subnetwork ii possesses unstable stationary states if the concentrations of the indicated species are sufficiently small. The reactions of ii are shown in thick lines in Figure 2. They provide a positive feedback loop regenerating autocatalytically NAD•, which, if coupled with the stable subnetworks, produces the oscillatory dynamcis. The instability generated by subnetwork ii is readily understood in terms of the topological features of its network diagram. At the steady state, the autocatalytic double-cycle through species coI, coII, and NAD• (reactions PO2, PO3, and PO4) has a zero β3 (a critical current cycle); there are no entrance reactions producing cycle species, and NAD• is removed by the exit reaction PO6, which is a topological arrangement that accounts for an instability.2 According to step ii in section 2, we can make a preliminary classification as follows: coI, coII, and NAD• are likely to be type X species (since they are on an autocatalytic cycle), O2 is type Y (due to exit step PO6), and β5 suggests that either per3+ or O2-, or both, are type W species. The species not indicated by any negative β can be either (essential) type Z or nonessential. Although helpful, this preliminary classification does not give a clue for the identification of type Z species and the amplitudes and shifts associated with oscillations must be examined in detail. Following our general scheme for calculations, we find the stationary state xs,j at the Hopf bifurcation and calculate the relative amplitudes of oscillations arel,j, the quench amplitudes qj, the phase shifts ∆φi,j of the jth species with respect to a selected ith species, and the concentration shifts sign(∂xs,i/∂x0,j) (sign symbolic representation is sufficient for categorization). The conservation constraint implies that the concentration shift
8560 J. Phys. Chem., Vol. 100, No. 20, 1996
Schreiber et al.
TABLE 2: Parameters of Model A and Reduced Model A at a Hopf Bifurcationc Section Aa
i per3+ coI coII coIII H2O2 NAD• O2O2
j
per3+
coI
coII
coIII
H2O2
NAD•
O2-
O2
arel,j qj ∆φNAD•,j sign(∂xs,i/∂x0,j)
0.2 16.0 -56
0.86 1 9
0.86 1.9 8
0.21 2.7 128
0.31 2.1 -23 + + + + + + -
1 2.1 0 + + + + + + -
0.63 3.0 19 + + + + + + -
0.3 7.7 162 + + + + + + +
Section Bb
i coI coII coIII H2O2 NAD• O2O2
j
coI
coII
coIII
H2O2
NAD•
O2-
O2
arel,j qj ∆φNAD•,j sign(∂xs,i/∂x0,j)
0.85 1 10 + + + + + -
0.85 2.0 10 + + + + + -
0.22 2.9 126 + + + + + -
0.32 2.2 -59 + + + + + -
1 2.0 0 + + + + + -
0.68 2.7 4 + + + + + -
0.34 7.4 161 + + + + + + +
The relative amplitudes arel,j, quench amplitudes qj, phase shifts ∆φNAD•,j with respect to NAD•, and sign symbolic concentration shifts sign(∂xs,i/ ∂x0,j) of model A at a Hopf bifurcation. For convenience, the relative amplitudes (quench amplitudes) are scaled so that the maximum (minimum) value equals 1. b The scaled relative amplitudes arel,j and quench amplitudes qj, the phase shifts ∆φNAD•,j, and sign symbolic concentration shifts sign(∂xs,i/∂x0,j) of a reduced model A (with fixed [per3+]) at a Hopf bifurcation. c Classification: coI, coII, NAD•stype X; O2stype Y; coIIIstype Z; H2O2snoness., O2-stype W. Category is 1CW. a
matrix is incomplete. The conditions imposed on the concentrations of the species coI, coII, NAD,˘ and O2 by the principal minor β4 and the consistency with available experimental data25,28 lead to the following choice of stationary concentrations: [per3+] ) 1 × 10-6 M, [coI] ) 0.01 × 10-6 M, [coII] ) 0.01 × 10-6 M, [coIII] ) 1 × 10-6 M, [H202] ) 0.1 × 10-6 M, [NAD•] ) 0.01 × 10-6 M, [O2-] ) 0.1 × 10-6 M, [O2] ) 1 × 10-6 M, [NADH] ) 20 × 10-6 M. The Hopf bifurcation occurs for any current Vs ) R1e1 + 3 j l ) Rl/∑l)1 Rl R2e2 + R3e3 such that the scaled coefficients R are R j ) 0.385, R j 2 ) 0.577, and R j 3 ) 0.038. Thus, the unstable current e2 contributes by 57.7% and dominates over the stable currents. The results are shown in Table 2A. NADH is nonessential by assumption, in accordance with experimental findings,18 and does not appear in the table. All other species have comparable relative amplitudes (within an order of magnitude) and hence might be labeled as essential. However, a closer examination is needed, since sometimes nonessential species, particularly if present in low concentrations, do not necessarily have small relative amplitudes and/or large quench amplitudes. We use two methods to identify such hidden nonessential species: (a) the Hopf bifurcation should not be significantly shifted if the steady state concentration value of a species suspected to be nonessential is fixed and (b) the Hopf bifurcation should be insensitive to an increase (even substantial) in the steady state concentration of that species (at fixed Rl’s). In effect, the first method amounts to setting the relative amplitude to zero and the quench amplitude to infinity. The second method makes the relative amplitude smaller and the quench amplitude larger as the concentration at steady state is increased. The quench amplitude of per3+ is significantly larger than the others, suggesting that it might be nonessential. This is corroborated by both methods. In addition, H2O2 is not part of the unstable extreme subnetwork ii, suggesting that its importance as a dynamical species for the occurrence of oscillations
is only marginal. Indeed, if the stationary value of [H2O2] is buffered or raised, the Hopf bifurcation is shifted negligibly, indicating that hydrogen peroxide is nonessential despite the “normal” values of the relative amplitude and the quench amplitude. Other species do not pass this test. Therefore, we shall assume that, in addition to NADH, per3+ and H2O2 are also nonessential.35 In the next step, the phase shifts are used to determine the role of essential species. The oscillations near the examined Hopf bifurcation display a nearly antiphase shift of oxygen with respect to NAD•, which itself is almost in-phase with coI and coII. This suggests that NAD•, coI, and coII are the autocatalytic species (type X) and O2 is the exit species (type Y). Compound III (coIII) advances NAD• but lags behind oxygen. Consequently, its role may be identified with the role of the feedback species (type Z). The O2- is only slightly in advance of NAD• and thus may be classified as either a type X or type W species. As pointed out in previous work,1,3 the methods of phase (or concentration) shifts cannot differentiate between type X and type W species. Based on the topology of the network in Figure 2, the role of O2- is that of the recovery species (type W) rather than the autocatalytic species, since O2- does not participate in the autocatalytic double-cycle. The large negative phase shift of per3+ is inconsistent with a possible shift behavior of any essential species and confirms its identification as a nonessential species. Notice that our preliminary classification based solely on the calculated negative β’s for the unstable current e2 is entirely consistent with that obtained from phase shifts. But the phase shifts calculations were necessary for identification of coIII as the feedback species. The second method of indicating the role of essential species is the analysis of the concentration shifts sign (∂xs,i/∂x0,j) in Table 2A. The columns (rows) corresponding to perturbations by (shifts of) the species NAD• and O2- are interchangeable, indicating the same or closely related type. Their self-regulation is normal suggesting that these are type X or type W species.
Oscillatory Enzymatic Reactions
J. Phys. Chem., Vol. 100, No. 20, 1996 8561
Compounds I and II (coI and coII) have the same shift behavior, whereas oxygen is inversely regulated with respect to all the previous species, which is consistent with type Y and with the 1C category. Compound III is negatively regulated with respect to all species but oxygen, and this consistent with type Z regulation for the 1C category even though the inverse selfregulation typical for type Z species is not observable because of the constraint. The shifts of coI and coII are consistent with type X for the 1C category. However, NAD•, O2-, coI, and coII are normally regulated with respect to O2, which differs from the prototypical shift behavior. A clearer picture emerges when per3+, indicated earlier are nonessential, is fixed. For the same steady state concentration values as in Table 2A, the Hopf bifurcation is found for a j 2 ) 0.594, and R j3 ) slightly shifted current with R j 1 ) 0.369, R 0.037. Table 2B shows that arel,j, qj, and ∆φNAD•,j are nearly the same as before except for the phase shift of H2O2, which is lagging NAD• significantly, pointing to the nonessential nature of hydrogen peroxide. Now the concentration shifts with respect to perturbations by all species are well defined, since the conservation constraint has been eliminated. The species coIII is inversely self-regulating, which confirms the suspected type Z behavior. Except for normal cross-regulation of type X and W species with respect to oxygen, the shift behavior coincides with the 1CW category prototype. In fact, we find that the 1CW prototype itself can be made to display either normal or inverse regulation of type X (or W) species with respect to type Y species. It is still possible to reduce the mechanism by fixing the concentration of the superoxide ion without losing the Hopf bifurcation. However, to see a Hopf bifurcation by fixing O2requires a current with R j 1 ) 0.141, R j 2 ) 0.845, and R j 3 ) 0.014, which differs significantly from the one for the nonreduced mechanism. Therefore, we conclude that model A belongs to the 1C category, switching between 1CX and 1CW subcategories depending on the external constraints. This is in agreement with the topological features of the network: there are entrance reactions producing type X species (reactions PO1 and PO7) characteristic of the 1CX category, but there is also an essential species (on the pathway formed by reactions PO6 and PO5) that behaves like a type X species but does not occur on the autocatalytic cycle. The mixing of the extreme subnetworks makes particular reactions more or less important, switching thereby the subcategory between 1CX and 1CW. A model considered by Fed’kina, Ataullakhanov, and Bronnikova27 is identical with a model A except for the termination reaction PO8, which is assumed to be of first order. This difference does not change the classification and categorization. 3.2. Model C. In an effort to find a model with complex oscillations seen in experiments, Aguda and Larter28 extended model A by adding a slow infusion of NADH and by adding a reaction pathway reducing per3+ by NAD• to a fifth form of the enzyme, per2+, which in turn is oxidized by O2 to coIII:
f NADH,
(PO11)
per3+ + NAD• f per2+ + NAD+,
(PO12)
per2+ + O2 f coIII.
(PO13)
The network diagram for the system of reactions PO1-PO13 is shown in Figure 3. The infusion makes [NADH] varying in time. The reduction of per3+ is reported to be slow for horseradish peroxidase,12 but it is faster when peroxidase from Coprinus cinereus is used.39 If PO12 and PO13 are considered to occur at a rate comparable to that of other reactions in the
(a)
(b)
Figure 3. Network diagrams of model C for the HRP reaction possessing two different unstable subnetworks.
mechanism,28 this additional loop seems to provide a period doubling route to chaos as observed in experiments.19,21 The network consists of 13 (pseudo)reactions PO1-PO13 and 10 dynamical species per3+, coI, coII, coIII, per2+, H2O2, NAD•, O2-, O2, and NADH. The enzyme species are bound by a conservation constraint, [per3+] + [coI] + [coII] + [coIII] + [per2+] ) const, and the steady state cone is four-dimensional. There are four extreme currents e1, ‚‚‚, e4 associated with the extreme subnetworks involving the following reactions: (i) PO9, PO10, (ii) PO2-PO6, PO9, PO11 (iii) PO1-PO3, PO6-PO9, PO11, and (iv) PO2-PO4, PO9, PO11-PO13. The subnetworks i, ii, and iii are identical with those for model A except for the additional infusion term PO11. Assuming mass action kinetics, the current e2 of subnetwork ii is unstable for precisely the same reasons as in model A. There is a second source of instability provided by the currents in the two-dimensional face generated by e3 and e4. The two unstable subnetworks are shown in thick lines in Figure 3. The difference between them is in the way per3+ is converted into coIII. The pathway PO12, PO13 makes it possible for per3+ to be first reduced to per2+ and the converted to coIII, rather than by a direct oxidation by O2-. Unlike with the first unstable subnetwork, the instability of the second one is no longer easy to explain in terms of topological features in its network diagram. Because of the stabilizing effect of reactions PO1 and PO8,2 there is no cycle with a zero β combined with an exit reaction, which is a typical source of instability.2 To gain insight, we examine a negative βk with the smallest k, which provides the simplest source of instability. We find that such a βk is dependent on the choice of current in the two-dimensional face formed by e3 and e4. In j 4e4, (where 0 < R j 3, R j4 < 1 particular, for all currents R j 3e3 + R
8562 J. Phys. Chem., Vol. 100, No. 20, 1996
Schreiber et al.
TABLE 3: Scaled Relative Amplitudes arel,j and Quench Amplitudes qj, Phase Shifts ∆ONAD•,j and Sign Symbolic Concentration Shifts sign (Dxs,i/Dx0,j) of Model C at a Hopf Bifurcation
i per3+ coI coII coIII per2+ H2O2 NAD• O2O2 NADH
j
per3+
coI
coII
coIII
per2+
H2O2
NAD•
O2-
O2
NADH
arel,j qj ∆φNAD•,j sign(∂xs,i/∂x0,j)
0.16 19.0 -52
0.82 1 2
0.82 1.9 1
0.18 2.7 132
1 3.3 -51
0.25 2.5 -30 + + + + + + -
0.95 2.1 0 + + + + + -
0.62 3.1 14 + + + + + -
0.29 8.4 162 + + + + + -
0.085 8400 98 + + + + + +
TABLE 4: Parameters of Model C and Reduced Model C at a Hopf Bifurcationc Section Aa j arel,j qj ∆φNAD•,j
coI
coII
coIII
per2+
H2O2
NAD•
O2-
O2
NADH
0.25 1 -34
0.25 2.7 -35
0.007 2.7 151
1 2.5 -27
0.22 1.6 27
0.21 1.6 0
0.31 1.8 87
0.82 37 146
0.15 370 91
3+
per
0.007 3.5 150
Section Bb i coI coII per2+ NAD• O2O2
j
coI
coII
per2+
NAD•
O2-
O2
∆φ sign(∂xs,i/∂x0,j) +
0 + + + + -
0 + + + + -
-22 + + + + -
0 + + + + -
99 + + + + -
154 + +
NAD•,j
a
The scaled relative amplitudes arel,j and quench amplitudes qj and phase shifts ∆φNAD•j of model C at a Hopf bifurcation. b The phase shifts ∆φNAD•,j and sign symbolic concentration shifts sign(∂xs,i/∂x0,j) of a reduced model C at a Hopf bifurcation. c Classification: coI, coII, NAD•stype X, O2stype Y, O2-stype Z, per2+stype W. Tentative category is 1CW.
and R3 + R4 ) 1), there is a negative β6 corresponding to the species per3+, per2+, H2O2, O2-, O2, and NADH, which must all be present in small concentrations to produce an unstable state. Such a condition is inconsistent with the available experimental data. In particular, NADH is typically present in j 4 is medium or large concentrations.18,22,26 However, if R sufficiently large, there is an additional negative β5 corresponding to the species coI, coII, NAD•, O2, and per2+. This condition is similar to that for the instability of subnetwork ii except that O2- is not among the species indicated by β5 (i.e., it might be of type Z or nonessential but not of type X, Y, or W) and per2+ is included. Thus, the autocatalytic double-cycle together with the exit reaction PO6 does support the instability, but in addition, [per2+] has to be small. This unstable feature involves the new pathway PO12, PO13 of model C, but because of the similarity with (ii), no drastic changes in the classification and categorization are expected. Nonetheless, the shifts have to be examined to draw definite conclusions about the classification. Provided that the unstable subnetwork ii (see Figure 3a) is dominant, the coupling of reactions PO12 and PO13 does not significantly alter conditioins for the occurrence of the Hopf bifurcation found for model A. We assume the same set of steady state concentrations as in section 3.1 and extend it with [per2+] ) 0.01 × 10-6 M. The current providing a Hopf j 2 ) 0.610, R j 3 ) 0.007, bifurcation is given by R j 1 ) 0.348, R and R j 4 ) 0.035 so that e2 is dominant. The results are listed in Table 3. The values of arel,i and qi indicate NADH and per3+ as nonessential. Both supplementary methods indicate that per2+ and H2O2 are nonessential as well. The phase shift relations of the essential species are similar to those in Table 2A, which suggests that the classification and categorization remain unchanged. By comparing data in Tables 2A and 3, we find
inconsistencies in the concentration shift matrix, which, however, entirely disappear upon reduction of the system by fixing [NADH] (not shown). This observation suggests caution in evaluating concentration shift data. A possible reason for the conflicting concentration shift behavior is that NADH and coIII can exchange their roles if NADH is present in small concentration and coIII is in excess. This competition alters the shift behavior unless one of the species is fixed. Next we consider a case of a dominant current in the unstable two-dimensional face associated with the subnetwork in Figure 3b. With the steady state values as above, a Hopf bifurcation j 2 ) 0.0009, R j 3 ) 0.0095, and R j4 ) occurs if R j 1 ) 0.0471, R 0.9425. As suggested by β5, we have chosen a current located near e4 in the two-dimensional formed by e3 and e4. Relative amplitudes, quench amplitudes, and phase shifts are shown in Table 4A. As in Table 3, the concentration shifts are inconsistent with those for the reduced models and therefore are omitted. The species per3+, coIII, and NADH are clearly indicated as nonessential. By fixing or raising its concentration, an additional nonessential species H2O2 is identified. The essential species coI, coII, and per2+ are in-phase, lagging NAD• by some 30°. This phase lag is small enough to justify an assumption that all these species are of the same or related type (presumably of type X or W). Oxygen is antiphase relative to coI, coII, and per2+, and the phase of O2- is in between. This is typical for the 1C category with a suggested classification of coI, coII, and NAD• as type X species, per2+ as a type W (it is off the autocatalytic double-cycle), oxygen as a type Y, and superoxide ion as a type Z species. The reduced model displays a Hopf bifurcation at a slightly shifted current with R j 1 ) 0.0397, j 3 ) 0.0079, and R j 4 ) 0.9516. The phase and R j 2 ) 0.0008, R concentration shifts of the reduced model shown in Table 4B
Oscillatory Enzymatic Reactions seem to confirm the suggested classification and categorization: the concentration shift matrix now coincides exactly with the prototypical shifts for the ICW category. None of the species in Table 4A can be fixed to obtain an even more reduced model displaying a Hopf bifurcation, albeit substantially shifted. However, the topological features of the network point to a more complex nature of the network’s instability. Namely, reaction PO7, through which O2- is expected to regulate autocatalytic production, does not provide a simple feedback Via tangent reaction PO4 as coIII does under conditions when it is essential. The negatiVe feedback comes about via the pathway PO7, PO12, PO13, PO6. Moreover, there are two unstable features within the network. The first one is the unstable feature contained in model A based on the autocatalytic double-cycle involving NAD•, coI, and coII destabilized by reaction PO6. The other involves NAD•, per2+, and O2 via reactions PO12, PO13, and PO6. The latter source of instability is not covered by any of the prototypes. There is no cycle, passing through any of the indicated species, that would be combined with (and destabilized by) an exit reaction. In this instability O2 is consumed by two competing reactions PO6 and PO13, and the two other reactants in these reactions, NAD• and per2+, are linked via reaction PO12. Such a feature provides a positive feedback via the pathway PO12, PO13, PO6, which is a basis for a competitive autocatalysis.38 The two unstable features are coupled and make this oscillator a mixture of two different instabilities (a nonsimple oscillator). Thus, model C offers two alternatives for classification and categorization. One is analogous to model A, that is, 1CW or 1CX category with coIII of type Z and O2- either of type W or nonessential, depending on the external constraints. The second one suggests a nonsimple oscillator with O2- as the feedback species. It may be viewed as a 1CW category oscillator with a complex negative feedback loop. Then coI, coII, and NAD• are type X species, per2+ is type W, and O2 is type Y. The second alternative assumes reaction PO12 to be sufficiently fast, which contradicts the earlier experimental data with horseradish peroxidase.36 Recently, experiments with peroxidase obtained from Coprinus indicate39 that per2+ is present in significant amounts. This seems to suggest that the pathway via PO12, PO13, and PO4 rather than via PO5 and PO4 might constitute the major negative feedback loop. However, our analysis shows that the network of model C does not provide such a possibility. If PO5 is made slow relative to PO12, and nonsimple oscillatory is found for which coIII is nonessential and the negative feedback loop does not involve PO12, PO13 followed by PO4. Rather, the competitive autocatalysis causes O2- to act as the feedback species via pathway PO7, PO12, PO13, PO6. Furthermore, in this case coIII and per3+ will oscillate in phase (see Table 4A), which contradicts the antiphase relation found experimentally both with horseradish10 and Coprinus39 peroxidases. We conclude that (a) the pathway via PO12, PO13, and PO4 cannot replace the pathway via PO5 and PO4 in regulating the autocatalytic loop and (b) the nonsimple oscillator is inconsistent with experimental data. 3.3. A Model Involving DCP. It is known experimentally that oscillations in the PO system occur in the presence of certain stabilizers, such as 2,4-dichlorophenol (DCP) and methylene blue (MB). Particularly, the first one is typically added in significant amounts to observe sustained oscillations experimentally, and therefore, it is desirable to include it in the reaction network. Unfortunately, there is not much experimental evidence on the chemistry of DCP in this reaction. Olsen37 conjectured that DCP reduces compound III and that the
J. Phys. Chem., Vol. 100, No. 20, 1996 8563
Figure 4. Network diagram of model A extended by reactions involving a DCP breakup.
generated radical DCP• produces an additional NAD•,
DCP + coIII f DCP• + coI,
(PO14)
(PO15) DCP• + NADH f DCP + NAD• For simplicity we combine these reactions with model A (see Figure 4). The network consists of 12 (pseudo)reactions PO1PO10, PO14, PO15 and 10 dynamical species, those of model A, DCP, and DCP•. There are two conservation constraints, [per3+] + [coI] + [coII] + [coIII] ) c1 and [DCP] + [DCP•] ) c2. The steady state cone is four-dimensional, formed by four extreme currents e1, ‚‚‚, e4. The first three are identical with those for model A, and the subnetwork iv associated with e4 involves reactions PO2, PO3, PO5, PO6, PO8, PO9, PO14, and PO15. The extreme current e2 is a basis for an oscillatory instability completely analogous to that for model A. With the same set of steady state concentrations extended by [DCP] ) 19 × 10-6 M and [DCP•] ) 1 × 10-6 M, the Hopf bifurcation j 2 ) 0.667, R j 3 ) 0.028, and R j 4 ) 0.027. occurs at R j 1 ) 0.278, R The relative amplitudes and quench amplitudes that DCP and DCP• are nonessential, and the additional species per3+ and H2O2 are indicated as nonessential by fixing or raising their concentrations. The shift data for the reduced model coincide with those in Table 28. This is expected, since subnetwork iv (the only subnetwork specific for the DCP model) is weakly coupled. Thus, the oscillator is categorized in the same way as that of model A. The additional DCP pathway brings an additional unstable subnetwork consisting of subnetworks iii and iv. Reactions PO14 and PO15 are participating in this subnetwork, and only reaction PO4, which is a part of the autocatalytic double-cycle, is not contained in it. By inspecting negative βk’s with the smallest k, we find a β4 associated with the species per3+, H2O2, O2-, and NAD•. It points once again to the competitive autocatalysis already found in model C. However, different species are involved in this case: per3+ is consumed by two competing reactions PO1, PO5 and the two other reactants in these reactions, H2O2 and O2-, are linked via reaction PO7. Since the small concentration requirement for per3+ (it is essential) is inconsistent with experimental data,10 we do not examine this instability in detail here. An alternative way for DCP breakdown is31 DCP + coIII f DCP• + per2+ + O2-
(PO14a)
Unstable features of the PO mechanism with reaction PO14a instead of PO14 remain the same as those described above.
8564 J. Phys. Chem., Vol. 100, No. 20, 1996
Schreiber et al. TABLE 5: Phase Shifts ∆ONAD•,j and Sign Symbolic Concentration Shifts sign(Dxs,i/Dx0,j) of a Reduced AD Model at a Hopf Bifurcationa i coIII NAD• O2
Figure 5. Network diagram of AD model for the HRP reaction.
Thus, the unstable network ii of model A is present. Competitive autocatalysis is present as well, since the feature which makes it dominant is the creation of NAD• via PO14 or PO14a. The latter is again inconsistent with experimental measurements, suggesting that per3+ is nonessential.10 Recently, efforts have been made to bring the mechanisms of the PO reaction closer to the results of experiments29-31 by including some reactions that have been known but omitted earlier for the sake of simplicity of the analysis. Among those are the reactions
2O2- + 2H+ f O2 + H2O2
(PO16)
NADH + MB+ f MBH + NAD+
(PO17)
MBH + O2 + H+ f MB+ + H2O2
(PO18)
NADH + H2O2 + H+ f 2H2O + NAD+ (PO19) whereas reaction PO7 has been left out.29,30 Reactions PO17, PO18 describe the catalytic effect of the other modifier, methylene blue. They can be merged if MB+ is nonessential. Without any detailed discussion we conclude that none of the reactions PO16-PO19 (nor the removal of PO7) bring about an entirely different unstable subnetwork. Thus, conclusions about the core oscillator based on the unstable subnetwork PO2-PO6, PO9 drawn repeatedly in this section hold true for the extended models as well. On the other hand, a good quantitative agreement with experimental results certainly requires inclusion of all relevant species and reaction steps even if they are classified as nonessential. 3.4. AD Model. The model by Alexandre and Dunford34 involves the following steps (Figure 5):
per3+ + O2 + NAD• f coIII + 2NAD•
(AD1)
coIII + NAD• f per3+ + 2NAD•
(AD2)
( ) f NAD•
(AD3)
NAD• f product
(AD4)
f O2, O2 f
(AD5) (AD6)
It may be viewed as a highly skeletonized and simplified model A with the concentration of NADH assumed to be constant and included into reaction rate coefficients. The species O2-, coI, and coII are left out of the network, and reactions around H2O2 are replaced by reaction AD3 forming NAD• from a constant internal source. Reaction AD1 resembles somewhat a combination of reactions PO5-PO7. Unlike reaction PO7 of model A, there are two NAD• produced by reaction AD1. Reaction AD2 is obtained by merging PO2, PO3, and PO4. The termination
j
coIII
NAD•
O2
∆φNAD•,j sign(∂xs,i/∂x0,j)
52 + + -
0 + -
109 + + +
a Classification: NAD•stype X, O stype Y, coIIIstype Z. Cat2 egory is 1CX.
reaction AD4 is first order only. Oxygen is formed from and removed to an external source as in model A. The network has six (pseudo) reactions AD1-AD6 and four dynamical species: per3+, coIII, NAD•, and O2. The enzyme species satisfy [per3+] + [coIII] ) const. The cone is threedimensional with three extreme currents e1, e2, and e3 that involve the following reactions: (i) AD5, AD6, (ii) AD3, AD4, and (iii) AD1, AD2, AD4, AD5. Assuming mass action kinetics, we find that there is no unstable subnetwork (neither extreme nor nonextreme) in the sense of a negative βk. This is not contradicting the SNA method of calculating β’s as indicators of an instability, since a negative β is a sufficient rather than a necessary condition2 (see also Appendix). The absence of any βk < 0 signifies that there can be no saddle stationary state and therefore no bistability (which contradicts experimental results32), but there still can be an unstable focus surrounded by a limit cycle as shown in ref 34. Our calculations indicate that subnetwork iii (thick lines in Figure 5), even though lacking a negative β, is the dominant one to achieve the oscillatory instability. The region of sustained oscillations is rather limited in the space of reaction rate coefficients,34 which is likely because of the condensed skeletal character of the model. Therefore, the stationary state concentrations we used with model A do not provide an instability in this case. Hence, we take the values corresponding to ref 34 to be [per3+] ) 4.751 × 10-6 M, [coIII] ) 1.249 × 10-6 M, [NAD•] ) 0.1369 × 10-6 M, and [O2] ) 6.576 × 10-6 M. A Hopf bifurcation j 2 ) 0.0004, R j 3 ) 0.7511. The occurs for R j 1 ) 0.2485, R calculated relative amplitudes, quench amplitudes, and phase shifts indicate that the only nonessential species is per3+, which oscillates in exactly the opposite phase with coIII because the conservation constraint involves just these two species. The phase and concentration shifts of the reduced model with per3+ fixed are given in Table 5. Clearly, NAD• is the only autocatalytic species coupling two different autocatalytic cycles. Compound III advances NAD• by about 50° and oxygen by about twice as much. Therefore, both coIII and O2 may be classified as feedback species, each one entering a different autocatalytic cycle. O2 is a suspected type Y species, since it is not well in-phase with coIII, the other suspected type Z species. Concentration shifts suggest that NAD• is of type X and O2 of type Y. But compound III lacks the inverse selfregulation to be classified as the missing type Z species. This ambiguity may be resolved by taking into account that NAD• does not have to be produced in stoichiometric excess by reaction AD1 to observe oscillations in the network. Assuming that one NAD• is produced rather than two by AD1, the phase shift of O2 is increased to about 130° so that it becomes more apparent that oxygen is actually a type Y species. If the stoichiometric coefficient is set to a value less than 1, then the phase shift is about 160° and the inverse self-regulation of coIII occurs. Thus, we conclude that this is a category 1CX oscillator with somewhat irregular shift behavior. To allow for bistability, observed in experiments,32 the model needs to be extended, for example, by splitting reaction AD1
Oscillatory Enzymatic Reactions
J. Phys. Chem., Vol. 100, No. 20, 1996 8565 TABLE 6: Phase Shifts ∆OX,j and Sign Symbolic Concentration Shifts sign(Dxs,i/Dx0,j) of a Reduced DOP Model at a Hopf Bifurcationa i X Y O2 a
Figure 6. Network diagram of DOP model for the HRP reaction.
into two, linked via an intermediate. This points again to its similarity to model A. 3.5. DOP Model. The last model of the PO system to be discussed was suggested at an early stage of studying the peroxidase-oxidase system by Degn, Olsen, and Perram32 (see also Olsen and Degn15,33) as a general system with a “quadratic branching” anticipated in the peroxidase-oxidase mechanism (see Figure 6):
NADH + O2 + X f 2X
(DOP1)
2X f 2Y
(DOP2)
NADH + O2 + Y f 2X
(DOP3)
X f product
(DOP4)
Y f product
(DOP5)
()fX f O2 O2 f f NADH
(DOP6) (DOP7) (DOP8) (DOP9)
In a variant of this mechanism33 O2 is left out of the reaction DOP1 and 3X are produced by DOP3 rather than 2X. Here X and Y are intermediates (possibly radicals) that have been associated with species differently by different authors: X ≡ NAD•, Y ≡ coIII,25 and X ≡ O2-, Y ≡ NAD•.40 The model reproduces qualitatively not only periodic oscillations and period doubling route to chaos44 but also, unlike the realistic models, mixed mode periodic and chaotic oscillations, which are observed experimentally.15 The network consists of nine reactions DOP1-DPO9 and four dynamical variables X, Y, O2, and NADH. No conservation constraint is present and the cone is five-dimensional with seven extreme currents. The corresponding subnetworks involve the reactions: (i) DOP7, DOP8; (ii) DOP4, DOP6; (iii) DOP2, DOP5, DOP6; (iv) DOP2, DOP3, DOP5, DOP7, DOP9; (v) DOP2-DOP4, DOP7, DOP9; (vi) DOP1, DOP4, DOP7, DOP9; (vii) DOP1, DOP2, DOP5, DOP7, DOP9. Network v is the only extreme subnetwork that provides unstable steady states (thick lines in Figure 6). There is just one negative principal minor, corresponding to X and Y. The species X and Y are linked by an autocatalytic cycle with negative β2 (a strong current cycle). SNA theory states that this instability is due to a combination of the second-order reaction DOP2 with the firstorder removal reaction DOP4.2 The unstable subnetwork of the DOP model is nearly identical to that of the Sel’kov model,41 which is essentially a prototype for category 2C oscillators1 (compare to Figure 1). Subnetwork v supports periodic oscillations when combined with stable subnetworks i-iv. Subnet-
j
X
Y
O2
∆φX,j sign(∂xs,i/∂x0,j)
0 + + -
-8 + + -
95 + + -
Classification: X, Ystype X, O2stype Z. Category is 2C.
work vi contains an additional loop, which is identified as a Lotka damped oscillator.24,42 Subnetwork vi alone does not produce unstable stationary states but allows for an oscillatory decay. A competition between (v) and (vi) leads to complex dynamics and provides an interpretation in chemical terms of quasiperiodicity and mixed mode oscillations24 observed in the model. The categorization has been discussed in prior work.3,4 Here, we use stationary state concentration consistent with those in section 3.1, [O2] ) 1 × 10-6 M and [NADH] ) 20 × 10-6 M, and we set [X] ) 0.01 × 10-6 M and [Y] ) 0.01 × 10-6 M. A Hopf bifurcation can be easily found, for example for R j1 ) j 3 ) 0.134, R j 4 ) 0.134, R j 5 ) 0.558, R j6 ) 0.134, R j 2 ) 0.014, R 0.013, and R j 7 ) 0.013. The relative amplitudes and quench amplitudes indicate that NADH is nonessential. For the same set of R j ’s, the reduced model obtained by fixing [NADH] displays the phase and concentration shifts as shown in Table 6. X and Y are mutually in-phase, advanced by O2. This suggests that X and Y are type X species and O2 a type Z species. The concentration shifts lead to the same conclusion. It is also possible to find a Hopf bifurcation if O2 is fixed instead of NADH, and then NADH plays the role of a type Z species. Thus, we conclude that the DOP model is a category 2C oscillator, as suggested by the unstable network, with X and Y as autocatalytic species and either oxygen or NADH as a feedback species, the other being nonessential depending on the constraints. The variant of the DOP model by Olsen33 is of the same category with the same role of species. Thus, the DOP model differs significantly from all other models in both the category and the role of oxygen. By comparison of the DOP and AD models, a similar double autocatalytic cycle is found but the enzyme cycle is absent in the DOP model. This comparison carries over to model A, suggesting that X in the DOP model plays a role similar to that of NAD• in the other models. However, as discussed in ref 34, it is virtually impossible to identify Y with any species in more realistic models. From the categorization standpoint, Y is of type X, which prompts a comparison with coI or coII of the autocatalytic double-cycle in model A. But if so, then oxygen should be left out of the reaction DOP3. Alternatively, it is possible to associate O2- with Y, since O2- is in-phase with NAD• in model A or C. But the cycle between NAD• and O2in those models does not produce an instability necessary for oscillations and cannot be identified with the cycle through X and Y in the DOP model. Hence, the categorization method rules out a consistent assignment of Y to some species occurring in the other models. The second-order reaction DOP2 is crucial for instability of the DOP model. However, it is possible to make the reaction DOP2 of first order and retain the oscillatory dynamics by inserting an intermediate to the Lotka subnetwork vi. This rearrangement also makes oxygen a type Y species and brings the DOP model closer to the other peroxidaseoxidase models. In concluding the section on the peroxidase-oxidase reaction, we suggest some experiments that help to determine the role of the species and thus to discriminate among reaction mech-
8566 J. Phys. Chem., Vol. 100, No. 20, 1996 anisms in which these species have different roles. The DOP model is different from the rest, largely because of the type Z role of oxygen (it is classified as type Y in all other models). A phase advancement or antiphase shift of O2 with respect to NAD• (the only consistently identifiable type X species in the DOP model) is one possibility, the negative vs positive selfregulating concentration shift of O2 is another one. The latter has been found to be positive in experiments,10 which makes the DOP model irrelevant despite its ability to reproduce qualitatively various types of experimentally observed dynamics.43,44 All types of experimentally found complex dynamics such as period doubling route to chaos19 and mixed mode oscillations20,21 are reproduced also by models described by some subset of steps PO1-PO19 (see refs 28 and 30 for period doubling and refs 31 and 45 for mixed mode oscillations). The AD model is consistent with model A but seems to be oversimplified. To discriminate between model A and model C, either per2+ should be identified as essential or nonessential or coIII should be identified as (essential) type Z or nonessential. Measured relative amplitudes indicate that coIII is essential of type Z10 and suggest that the additional pathway in model C is not crucial. The determination of the role of the modifiers DCP and MB rests on whether they are essential or not. The experimental oscillations are observed traditionally in the presence of both DCP and MB but also in the absence of DCP29 or MB.39 In the models, both modifiers are nonessential for oscillations. This can be verified experimentally by measurements of relative amplitudes and/or quench amplitudes. Reactions PO14, PO14a, PO15 and PO17, PO18 serve the same purpose in the mechanism: to create initially some NAD•, which is an essential ingredient for oscillations to start (NAD• is type X species). This is accomplished directly according to PO14 or PO14a, or indirectly via PO17, PO18 followed by PO1. Thus, according to model mechanisms, the presence of modifiers should have a significant impact on the transient behavior before oscillations set in, which can be verified experimentally. 4. Conclusions The analysis given here for each proposed reaction mechanism of an oscillatory enzymatic reaction leads to extensive insight. We identify essential and nonessential species, the role of each essential and nonessential species, the connectivity of these essential species (the category), and the particular elementary reaction steps leading to an oscillatory instability. The study proceeds from proposed reaction mechanisms to predictions of various experiments such as relative amplitudes, quench amplitudes, phase shifts, and sign symbolic concentration shifts. Wherever possible we use these predictions to compare with experiments in order to discriminate among competing models of reaction mechanisms. The reverse procedure, that of deducing details of reaction mechanisms from these designed experiments, has been demonstrated for the inorganic chloriteiodide reaction4,9 and the horseradish peroxidase-oxidase reaction.10,31 Acknowledgment. This work has been supported in part by the National Science Foundation and the Department of Energy/Basic Energy Sciences. Supporting Information Available: Details of the analysis for reactions 2-5 and of the methods determining the categories
Schreiber et al. (50 pages). Ordering information is given on any current masthead page. References and Notes (1) Eiswirth, M.; Freund, A.; Ross, J. AdV. Chem. Phys. 1991, 90, 127; J. Phys. Chem. 1991, 95, 1294. (2) Clarke, B. L. AdV. Chem. Phys. 1980, 43, 1. (3) Chevalier, T.; Schreiber, I.; Ross, J. J. Phys. Chem. 1993, 97, 6776. (4) Stemwedel, J. D.; Schreiber, I.; Ross, J. AdV. Chem. Phys. 1995, 89, 327. (5) Strasser, P.; Stemwedel, J. D.; Ross, J. J. Phys. Chem. 1993, 97, 2851. (6) Finkeova, J.; Dolnik, M.; Hrudka, B.; Marek, M. J. Phys. Chem. 1990, 94, 4110. (7) Hynne, F.; Sørensen, P. G.; Nielsen, K. J. Chem. Phys. 1990, 92, 1747; J. Chem. Phys. 1990, 92, 4778. (8) Vance, W.; Ross, J. J. Chem. Phys., accepted for publication. (9) Stemwedel, J.; Ross, J. J. Phys. Chem. 1995, 99, 1988. (10) Hung, Y.-F.; Ross, J. J. Phys. Chem. 1995, 99, 1974. (11) Arkin, A.; Ross, J. J. Phys. Chem. 1995, 99, 970. (12) Yamazaki, I.; Yokota, K.; Nakajima, R. Biochem. Biophys. Res. Commun. 1965, 21, 582. (13) Yamazaki, I.; Yokota, K. Biochim. Biophys. Acta 1967, 132, 310. (14) Degn, H. Biochim. Biophys. Acta 1969, 180, 271. (15) Olsen, L. F.; Degn, H. Nature 1977, 267, 177. (16) Aguda, B. D.; Hofmann Frisch, L.-L.; Olsen, L. F. J. Am. Chem. Soc. 1990, 112, 6652. (17) Lazar, J.; Ross, J. Science 1990, 247, 189. (18) Samples, M.; Hung, Y.-F.; Ross, J. J. Phys. Chem. 1992, 96, 7338. (19) Geest, T.; Steinmetz, C. G.; Larter, R.; Olsen, L. F. J. Phys. Chem. 1992, 96, 5678. (20) Hauck, T.; Schneider, F. W. J. Phys. Chem. 1993, 97, 391. (21) Hauck, T.; Schneider, F. W. J. Phys. Chem. 1994, 98, 2072. (22) Yokota, K.; Yamazaki, I. Biochemistry 1977, 16, 1913. (23) Aguda, B. D.; Clarke, B. L. J. Chem. Phys. 1987, 87 (6), 3461. (24) Aguda, B. D.; Larter, R.; Clarke, B. L. J. Chem. Phys. 1989, 90, 4168. (25) Aguda, B.; Larter, R. J. Am. Chem. Soc. 1990, 112, 2167. (26) Olsen, L. F.; Degn, H. Biochem. Biophys. Acta 1978, 523, 321. (27) Fed’kina, V. R.; Ataullakhanov, F. I.; Bronnikova, T. V. Biophys. Chem. 1984, 19, 259. (28) Aguda, B.; Larter, R. J. Am. Chem. Soc. 1991, 113, 7913. (29) Olson, D. L.; Williksen, E. P.; Scheeline, A. J. Am. Chem. Soc. 1995, 117, 2. (30) Bronnikova, T. V.; Fed’kina, V. R.; Schaffer, W. M.; Olsen, L. F. J. Phys. Chem. 1995, 99, 9309. (31) Hung, Y.-F.; Schreiber, I.; Ross, J. J. Phys. Chem. 1995, 99, 1980. (32) Degn, H.; Olsen, L. F.; Perram, J. Ann. N. Y. Acad. Sci. 1979, 316, 623. (33) Olsen, L. F. Phys. Lett 1983, 94A, 454. (34) Alexandre, S.; Dunford, H. B. Biophys. Chem. 1991, 40, 189. (35) Olson, D. L.; Scheeline, A. J. Phys. Chem. 1995, 99, 1204; J. Phys. Chem. 1995, 99, 1212. (36) Wittenberg, J. B.; Noble, R. W.; Wittenberg, B. A.; Antonini, E.; Brunori, M.; Wyman, J. J. Biol. Chem. 1967, 242, 626. (37) Olsen, L. F. Personal communication. (38) Eiswirth, M. Suri Kagaku 1994, 372, 59. (39) Kummer, U.; Valeur, K. R.; Baier, G.; Wegmann, K.; Olsen, L. F. Biochim. Biophys. Acta, submitted. Valeur, K. R.; Olsen, L. F. Biochim. Biophys. Acta, submitted. (40) Olsen, L. F.; Hofman Frisch, L.-L. H.; Schaffer, W. M. In A. Chaotic Hierarchy; Baier, G., Klein, M., Eds.; World Scientific Publishers: Singapore, 1991; p 299. (41) Selkov, E. E. Eur. J. Biochem. 1968, 4, 79. (42) Lotka, A. J. Phys. Chem. 1910, 14, 271. Volterra, V. Mem. Acad. Lincei 1926, 2, 31. (43) Steinmetz, C. G.; Larter, R. J. Chem. Phys. 1991, 94, 1388. (44) Steinmetz, C. G.; Geest, T.; Larter, R. J. Phys. Chem. 1993, 97, 5649. (45) Fed’kina, V. R.; Bronnikova, T. V.; Ataullakhanov, F. I. Biophysics (Engl. Transl.) 1992, 37, 781.
JP952853X
