J. Phys. Chem. 1996, 100, 18911-18915
18911
Reduction of Dimension of a Chemically Realistic Model for the Peroxidase-Oxidase Oscillator Daniel S. Bensen and Alexander Scheeline* Center for Complex Systems Research, Beckman Institute for AdVanced Science and Technology, UniVersity of Illinois at Urbana Champaign, 405 North Mathews, Urbana, Illinois 61801 ReceiVed: December 6, 1995; In Final Form: March 11, 1996X
We develop a low-dimensional model of the peroxidase-oxidase oscillator based on the higher-dimensional Urbanalator model. The low-dimensional model is consistent with the chemical foundations of the Urbanalator, accurately reproduces its oscillatory behavior, and also bears some similarities to the simpler abstract models. The remaining differences may suggest directions for further model refinements.
1. Introduction
TABLE 1: Reactions in the Urbanalator
The peroxidase-oxidase reaction has been an object of research for many years. Its dynamical richness was first brought to light by Yamazaki et al.1 In parallel with experiments characterizing chemical oscillations in this reaction among O2, NADH, a peroxidase enzyme typically from horseradish root but on occasion from bovine milk2 or other sources,3 and modifiers such as 2,4-dichlorophenol and methylene blue, there have been significant theoretical efforts to model the chemistry involved. Some of these models have been abstract,4,5 involving four species A, B, X, and Y whose connections to real chemicals were at times tenuous. Others6-11 have been chemically specific. Until quite recently, the correspondence between chemical models and observed dynamics was inadequate for predictive purposes. The ABXY models, in contrast, while illsuited to quantitative prediction, gave excellent qualitative agreement with observed oscillatory behavior. It thus seemed plausible that some limiting case of a chemically realistic model should recover the dynamics described in the ABXY models, using fewer variables than the nine9 or ten10,13 in the chemically realistic cases. Toward that end, we have started with our nine-species model and simplified it by noting intermediate steps that can be combined into approximate unified steps. This analysis gives additional insight into the functioning of the chemical models and features of those models which allow for oscillation to occur. The oscillation of the peroxidase-oxidase system is driven by oxidation of NADH by aqueous oxygen. It is characterized by an accumulation of Per6+ and aqueous oxygen, followed by consumption of oxygen, and finally by conversion of the Per6+ back to Per3+. The last step occurs almost instantaneously. After the conversion, there is a short refractory period during which excess NAD• dimerizes (see Figure 2). While Per2+ is absent from the current version of the Urbanalator, it has been incorporated in a recent model derived from the Urbanalator10 and may be used in the future. It may play a role similar to those of Per4+ and Per5+, which are eliminated in the final low-dimensional model, but determine the model dynamics through the chemically allowable stoichiometry. Several species are conspicuously absent during all or most of the oscillation: Per4+, Per5+, NAD•, O2-•, and H2O2. Also, the concentration of Per3+ is relatively low late in the oscillation. Each of these species is treated as being in a quasi steady state X
Abstract published in AdVance ACS Abstracts, June 1, 1996.
S0022-3654(95)03611-2 CCC: $12.00
R1: R2: R3a: R3b: R4: R5a: R5b: R6: R7: R8: R9a: R9b: R10:
NADH + O2 + H+ f NAD+ + H2O2 Per3+ + H2O2 + 2H+ f Per5+ + 2H2O Per5+ + NADH f Per4+ + NAD• + H+ Per4+ + NADH f Per3+ + NAD• + H+ NAD• + O2 f NAD+ + O2-• Per3+ + O2-• + 4H+ f Per6+ + 2H2O 2O2-• + 2H+ f H2O2 + O2 Per6+ + NAD• f Per5+ + NAD+ NADHR f NADH NADH f NADHp O2(g) f O2(aq) O2(aq) f O2(g) 2NAD• f (NAD)2
(QSS). This allows us to set the time-evolution function for each species to zero, resulting in an equation that places a constraint on the reaction rates involved. Each constraint is a stoichiometric relationship which can be used to construct a set of net reactions simpler than the original system. The constraining equations also provide the actual QSS concentrations. In the case of O2-•, there are two distinct states, or phases, which are treated separately. What remains in an interaction between NADH, aqueous oxygen, hydrogen peroxide, and one oxidation state of peroxidase, for which either Per3+ or Per6+ can be used. 2. Notation and Definitions 2.1. Notation Used in This Paper. The Urbanalator is defined by a set of chemical reactions which determine the timeevolution functions of the involved species. The reactions will be written as Rn, where n is the number or label of the reaction; rn will denote the rate of Rn. In the time-evolution equations, the time rate of change dX/dt of any concentration X will be written as X˙ , and the species concentrations will be abbreviated as follows: A ≡ [O2(aq)], B ≡ [NADH], H ≡ [H2O2], N ≡ [NAD•], O ≡ [O2-•], and Pn ≡ [Pern+]. n is the peroxidase oxidation state and may equal 3, 4, 5, or 6 in this model. Where two or more reactions are grouped together in a series to define a net reaction, the reaction symbols will be enclosed between braces, as for instance R12 ≡ {R1,R2}. In the equation for the series reaction, its rate will be written above the arrow in the equation. 2.2. Definition and Behavior of the Urbanalator. The reactions in the full Urbanalator are given in Table 1. NADHR is a reservoir of NADH, and NADHp refers to unreactive degradation products (for example, anomerization may lead to unreactive species). The rate of each reaction is knX1nX2n, where © 1996 American Chemical Society
18912 J. Phys. Chem., Vol. 100, No. 49, 1996
Bensen and Scheeline
Figure 1. Concentrations during one cycle. (a, left) Non-QSS species throughout the oscillation. B is scaled down by a factor of 10; H can barely be seen starting at t ) 400 s. (b, middle) H and O increase in phase II. (c, right) Per4+, Per5+, and NAD• appear during the transition.
TABLE 2: Olsen’s 1983 ABXY Model R1: R2: R3: R4:
B + X f 2X 2X f 2Y A + B + Y f 3X X f products
R5: R6: R7: R8:
Y f products X0 f X A0 T A B0 f B
kn is the rate constant for Rn, and X1n and X2n, if present, are the reactant concentrations, excluding H+. For a discussion of the numerical values used, see Olson et al.9 In R7 and R9a, the (constant) reactant concentrations are absorbed into the rate constants. Ignoring the unreactive species (NAD)2, the time-evolution equations are
A˙ ) -k1AB - k4AN + k5bO + k9a - k9bA
(2.1)
B˙ ) -k1AB - k3aP5B - k3bP4B + k7 - k8B
(2.2)
H˙ ) k1AB - k2P3H - k5bO2
(2.3)
2
species with the same letter, but in italics. The time-evolution equations are then
A˙ ) k7(A0 - A) - k3ABY
(2.10)
B˙ ) k8B0 - k1BX - k3ABY
(2.11)
X˙ ) k1BX - 2k2X2 + 3k3ABY - k4X + k6X0 (2.12) Y˙ ) 2k2X2 - k3ABY - k5Y
(2.13)
In comparing the Urbanalator with Olsen’s model, the nth reaction in the former will be written as R(U) n , that in the latter as R(O) n . 3. Simplification of the Urbanalator
N˙ ) k3aP5B + k3bP4B - k4AN - k6P6N - 2k10N
2
(2.4)
O˙ ) k4AN - k5aP3O - 2k5bO2
(2.5)
P˙ 3 ) -k2P3H + k3bP4B - k5bO2
(2.6)
P˙ 4 ) k3aP5B - k3bP4B
(2.7)
P˙ 5 ) k2P3H - k3aP5B + k6P6N
(2.8)
P˙ 6 ) k5aP3O - k6P6N
(2.9)
3.1. Preliminary Simplification. These reaction pairs always occur together in the model: R7 and R8, and R9a and R9b. Therefore, we will combine each pair into a pseudoreaction which contains the dynamical properties of both reactions.
RB: NADHother f NADH RA: O2(g) f O2(aq)
Figures 1 and 2 show the concentrations and reaction rates during one oscillation of the Urbanalator. Stoichiometric relations can be determined from the reaction rates by noting which curves overlap or are otherwise correlated. 2.3. Definition of One Abstract Model. Among the abstract ABXY models, Olsen’s 1983 model5 describes the known PO dynamics accurately and concisely. We will use this model for comparison with the reduced version of the Urbanalator. The reactions in Olsen’s model are given in Table 2. k1 is associated with peroxidase, A with oxygen, B with NADH, X usually with NAD•, and Y possibly with O2-• or Per6+. Following Olsen, we represent the concentration of each
The rates of these reactions are rB ) r7 - r8 and rA ) r9a - r9b. The “other” subscript in RB refers to both reservoir NADH and NADH products, since we are only concerned with the concentration of NADH. (NAD)2 is assumed to be unreactive and is therefore ignored. 3.2. Elimination of Per4+ and Per5+. Per4+ and Per5+ are absent during virtually the entire oscillation, implying that they are used up as soon as they are produced. This fact, which is due to the high concentration of NADH and the moderate values of the reaction constants, means that the reactions that consume these species are at no time rate-limiting. Therefore, we may incorporate Per4+ and Per5+ in combined series reactions involving the two reactions that produce Per5+: R2 and R6. By combining each of these reactions with R3a and R3b, we will compress four reactions into two, and eliminate P4 and P5 as dynamical variables. Defining RL ≡ {R2,R3a,R3b}, we find
Peroxidase-Oxidase Oscillator Model
J. Phys. Chem., Vol. 100, No. 49, 1996 18913
Figure 2. Reaction rates during one cycle. The change in stoichiometry is evident midway through the oscillation, the rates increase dramatically before and during the transition, and the refractory period can be seen as a shoulder immediately after.
that RL, in which peroxidase loops through its three lowest oxidation states, is an engine for NAD• production.
R2: Per3+ + H2O2 + 2H+ f Per5+ + 2H2O R3a: Per5+ + NADH f Per4+ + NAD• + H+ R3b: Per4+ + NADH f Per3+ + NAD• + H+ r
RL: 2NADH + H2O2 928 2NAD• + 2H2O Note that Per3+ is eliminated from the stoichiometry, although it affects the kinetics. The path by which Per6+ decays into Per3+ is R63 ≡ {R6,R3a,R3b}.
R6: Per6+ + NAD• f Per5+ + NAD+ R3a: Per5+ + NADH f Per4+ + NAD• + H+ R3b: Per4+ + NADH f Per3+ + NAD• + H+
The only oxygen-consuming reaction in Olsen’s model is (U) (U) (U) R(O) 3 , which is similar to R1L ≡ {R1 ,RL }, 3NADH + O2 + H+ f 2NAD• + NAD+ + 2H2O. However, R(O) 3 dominates during small-k1 simulations, which correspond to low peroxidase concentration; it also turns off when oxygen (A) accumulates, (O) indicating peroxidase activity through R(O) 1 . Therefore, R3 seems to represent reactions, at least some of which are independent of peroxidase. For an alternate interpretation of (U) (U) R(O) 3 , R5a could be combined somewhat arbitrarily with R1L to produce the net reaction O2 + 3NADH + Per3+ + O2-• f Per6+ + 2NAD•. Two notable differences between the two models are that Olsen includes degradation of Y by a first-order reaction to (O) unreactive products (R(O) 5 ) and an external source of X (R6 ). -• In constrast, the Urbanalator gives degradation of O2 by a second-order reaction to give reactive product H2O2 and includes no external source of NAD•.
r
R63: Per6+ + 2NADH 968 Per3+ + NAD• + NAD+ + 2H+ R2 and R6 are assumed to be rate-limiting, so that rL ) r2 and r63 ) r6. This means that r3a ) r3b ) r2 + r6. The second equals sign gives P4 ) (k2P3H + k6P6N)/k3bB, and the first yields P5 ) (k3b/k3a)P4. These are the QSS concentrations of Per4+ and Per5+. 3.3. Parallels with Olsen’s 1983 Model. The NADH (B) feed in this model is very similar to that of the Urbanalator, and the oxygen (A) feeds of the two models are identical in form. (O) (U) (U) As a group, R(O) 2 and R4 resemble R4 and R10 : Each pair • has one reaction that converts X(NAD ) to Y(possibly O2-•) and one reaction quadratic in X(NAD•). The terms are not combined in the same way in the two models but may produce similar behavior nonetheless. The clearest correlation among the core oscillatory reactions (O) is between the R(U) 63 and R1 , confirming the fairly wellaccepted association of X with NAD•. Not only are both reactions autocatalytic in NAD• (X) but both consume NADH (B) and both depend on peroxidase (k1).
4. Stoichiometry of the Urbanalator The reactions in the Urbanalator can be further combined to show their stoichiometric relationships in periodic oscillations. When a reaction is removed from the system, it is replaced by a switch which changes either the reactant concentrations or the stoichiometry. This is done whenever an approximation made in combining the reactions becomes invalid. 4.1. Elimination of NAD•. NAD• is also absent through most of the reaction cycle, due to the large reaction constant of R4, which consumes NAD•. R4 can be absorbed into series reactions, since it is never rate-limiting, eliminating NAD• as an explicit dynamical variable. NAD• quickly reaches a QSS concentration, which is determined from the NAD• timeevolution equation:
N˙ ) 2k2HP3 - (k4AN - k6P6N) ) 2k2HP3 - DN
(4.1) (4.2)
with D ≡ k4A - k6P6. N˙ ) 0 gives N ) 2k2HP3/D. RL, R63,
18914 J. Phys. Chem., Vol. 100, No. 49, 1996 and R4 are compressed into two series reactions, RL4 ≡ {RL, 2R4} and R64 ≡ {R63,R4}.
RL:
r2
2NADH + H2O2 98 2NAD• + 2H2O
2R4: 2NAD• + 2O2 f 2NAD+ + 2O2-• r
RL4: 2NADH + H2O2 + 2O2 928 2NAD+ + 2O2-• + 2H2O r
R63: Per6+ + 2NADH 968 R4:
Per3+ + NAD• + NAD+ + 2H+ NAD + O2 f NAD+ + O2-•
Bensen and Scheeline These two stoichiometries continue until Per3+ is consumed, at which point H2O2 and R5b both appear. 5.2. Phase II: Consumption of O2(aq). Whereas R1 limited the overall reaction rate in phase I, the concentrations of O2 and H2O2 determine the reaction rates in phase II. H2O2 accumulates (see Figure 1b), keeping Per3+ in a small QSS concentration. Setting P˙ 3 ) 0 yields r64 ) r36, implying that L1 continues in phase II. In addition to L1, the QSS condition for O2-• requires that 2rL4 + r64 ) r36 + 2r5b, or rL4 ) r5b. See Figure 2 for the change in stoichiometry. The time-evolution equations for O and P3 are O˙ ) r4 - r36 - 2r5b and P˙ 3 ) r63 - r36, or
•
r6
O˙ ) (k4A)N - k5aOP3 - 2k5bO2
(5.1)
P˙ 3 ) (k6P6)N - k5aOP3
(5.2)
2k6P6k2H P3 ) k5aOP3 D
(5.3)
R64: Per6+ + 2NADH + O2 98 Per3+ + 2NAD+ + 2H+ + O2-•
P˙ 3 ) 0 leads to
The net effect of RL4 and R64 is to oxidize NADH, producing O2-•; R64 also converts Per6+ into Per3+. 5. Elimination of O2-•. Phases of the Oscillation. Like NAD•, O2-• is a free radical and reacts rapidly; therefore, it may also be eliminated as an explicit model variable. However, there are now two net reactions that produce O2-•, and so this process will not simplify the model as much as was the case for NAD•. We can, however, achieve good simplification by treating the accumulation and depletion phases separately. 5.1. Phase I: Accumulation of Per6+ and O2(aq). Phase I is characterized by small QSS concentrations of hydrogen peroxide and superoxide radical, allowing us to set their timeevolution functions to zero. In addition, due to the small concentration of superoxide, and moderate quantities of Per3+, R5a overwhelms R5b in the competition for O2-•. This renders R5b negligible and allows us to ignore R5b during phase I. Using the QSS condition for H2O2 and the absence of R5b, we find that rL ) r1; in other words, R1 limits the rate of RL4. Then k2HP3 ) r1, or H ) r1/(k2P3). Also, N ) 2r1/D during phase I. The stoichiometry consistent with this condition is a net producer of Per6+, N1 ≡ {R1, RL4, 2R36}.
R1:
NADH + O2 + H+ f NAD+ + H2O r
RL4: 2NADH + H2O2 + 2O2 928 2NAD+ + 2O2-• + 2H2O 2R5a: 2Per3+ + 2O2-• + 8H+ f 2Per6+ + 4H2O N1:
r1
3NADH + 3O2 + 2Per3+ + 9H+ 98 3NAD+ + 2Per6+ + 6H2O
The rate of N1 is rN1 ) r1 ) k1AB. Doing the same thing for superoxide radical gives us r4 ) r36, meaning that Per3+ is converted to Per6+ by any O2-• produced by R4. k4AN ) k5aP3O gives us O ) k4AN/(k5aP3) ) 2k4Ar1/(k5aP3D). The required stoichiometry is L1 ≡ {R64, R36}. r
R64: Per6+ + 2NADH + O2 968 3+
R36: Per
+ O2
-•
Per3+ + 2NAD+ + 2H+ + O2-• + 4H+ f Per6+ + 2H2O r
L1: 2NADH + O2 + 2H+ 968 2NAD+ + 2H2O The rate of L1 is rL1 ≡ r6 ) k6P6N ) (2k6P6r1)/D.
which, after eliminating P3, gives O ) cOH/D with cO ≡ 2k2k6P6/ k5a. Finally, setting O˙ ) 0 yields
2k5bO2 )
[
]
2k4Ak2H 2k2k6P6H P3 D D
) 2k2HP3
(5.4) (5.5)
Then we have P3 ) k5bO2/(k2H) and N ) 2k5bO2/D. The O2-• QSS requires a new series reaction, L2 ≡ {R5b, RL4}. It is
R5b: 2O2-• + 2H+ f H2O2 + O2 RL4: 2NADH + H2O2 + 2O2 f 2NAD+ + 2O2-• + 2H2O L 2:
r5b
2NADH + O2 + 2H+ 98 2NAD+ + 2H2O
Note that this reaction has the same stoichiometry as L1 but has a different reaction rate, rL2 ) r5b ) k5bO2. Also, neither phase II stoichiometry is a net producer or consumer of H2O2, leaving R1 to accumulate H2O2 independent of the two loop reactions. The rates of L1 and L2 increase throughout phase II for two reasons. Accumulation of H2O2, which is in the numerator of the QSS concentrations, dominates initially. Once consumption of oxygen exceeds its supply, the denominator D starts to shrink. At the end it disappears precipitously, triggering the autocatalytic avalanche of NAD• via R63. 5.3. Transition and Refractory Period. As phase II proceeds, oxygen consumption and H2O2 concentration increase until D decreases to zero. At this point R4 shuts off, allowing NAD• to build up immediately as R64 suddenly converts all of the Per6+ back to Per3+ (see Figure 1c). This essentially instantaneous transition marks the end of the peroxidase cycle and begins a short refractory period which lasts approximately 1 s. During this period, R10 dominates due to the unusually large concentration of NAD• (see Figure 2). R2 shuts off, since NAD• also scavenges oxygen, and R4 and R36 continue at a fairly uniform rate, limited by production of oxygen by RA. 5.4. Low-Dimensional Model. To summarize, the final reduced model begins in phase I with RA supplying NADH, RB building up oxygen, L1 retarding that buildup somewhat, and N1 converting Per3+ into Per6+. The non-QSS species are oxygen, NADH, and peroxidase. The time-evolution equations
Peroxidase-Oxidase Oscillator Model
J. Phys. Chem., Vol. 100, No. 49, 1996 18915 abundant species. These reduced models shed light on connections between PO chemistry and the abstract ABXY models and clarify the stoichiometries involved in the chemical switch dynamics of the PO oscillator. In the reduction process, series reactions were introduced which may be used in further refinements of the Urbanalator. Acknowledgment. This work was supported in part by NSF Grant CHE-93-07549. Appendix
Figure 3. Oxygen concentrations during one cycle. A rises and falls smoothly during one oscillation, making it a good tool for comparing models. The reduced model runs faster than the Urbanalator during phase II.
of their concentrations are A˙ ) rA - 3rN1 - rL1, B˙ ) rB - 3rN1 - 2rL1, and P˙ 6 ) 2rN1, or
A˙ ) rA - k1AB(3 + (2k6P6/D))
(5.6)
B˙ ) rB - k1AB (3 + (4k6P6/D))
(5.7)
P˙ 6 ) 2r1
(5.8)
After complete oxidation of Per3+, RA and RB continue in phase II, while the built-up oxygen is completely consumed via L1 and L2. The non-QSS species are oxygen, NADH, and H2O2, and the time-evolution equations of their concentrations are A˙ ) rA - r1 - rL1 - rL2, B˙ ) rB - r1 - 2rL1 - 2rL2, and H˙ ) r1:
A˙ ) rA - k1AB - k6Ptot - k5b(cOH/D)2
(5.9)
B˙ ) rB - k1AB - 2k6Ptot - 2k5b(cOH/D)2
(5.10)
H˙ ) k1AB
(5.11)
The oscillation is then ended by the essentially instantaneous reduction of Per6+ by R63, followed by a very brief interlude dominated by R10. As Figure 3 shows, the phase I dynamics is essentially exact; phase II has the correct qualitative shape, although the reduced model runs faster than the full model. To complete the model, a discreet variable, which we will name φ, is introduced to indicate the phase. φ may arbitrarily be set to equal 0 during phase I and 1 during phase II. The entire model is then four-dimensional: The first variable is A; the second variable is B; the third, which we will call C, represents P6 during phase I and H during phase II; and the fourth is φ. When φ ) 0 and C reaches Ptot, φ is set to 1, C is set to 0, and eqs 5.6, 5.7, and 5.8 are used. When φ ) 1 and D ≡ k4A - k6Ptot decreases to 0, B is decreased by an amount 2(Ptot + A + C); A, C, and φ are set to 0; and eqs 5.9, 5.10, and 5.11 are used. 6. Conclusion The Urbanalator can be approximated by lower-dimensional models derived from QSS approximations for reactive and/or
Stability of the QSS Concentrations. Within the QSS approximation (assuming the QSS concentrations are dynamical fixed points), the stability analysis for phase I is straightforward. We define δH ≡ H - Hqss, δN ≡ N - Nqss, and δO ≡ O - Oqss, where the “qss” subscripts refer to QSS values. After linearizing around the QSS point for H, we get
δ˙ H ) -(k2P3)δH
(A.1)
which is just the equation for a decaying exponential with time constant k2P3. This means that δH f 0. The linearized equation for N is
δ˙ N ) 2(k2P3)δH - DδN
(A.2)
Since δH f 0, this equation approaches that of a decaying exponential (time constant D), and so δN f 0. Finally, the same thing happens for O:
δ˙ O ) (k4A)δN - (k5aP3)δO
(A.3)
implies δO f 0, time constant k5aP3. Thus the QSS point is stable. In phase II, the stability analysis is more difficult. The equations are more coupled, resulting in a cubic equation with coefficients which change in time. This equation must therefore be solved many times throughout phase II and will not be very accurate toward the end of the oscillation, when the QSS approximation breaks down. O2-• does in fact break out of QSS part way through phase II, slowing down the overall reaction. The oscillation as a whole, however, is fairly accurate. References and Notes (1) Yamazaki, I.; Yokota, K.; Nakajima, R. Biochem. Biophys. Res. Commun. 1965, 21, 582-586. (2) Nakamura, S.; Yokota, K.; Yamazaki, I. Nature 1969, 222, 794. (3) Kummer, U.; Valeur, K. R.; Baier, G.; Wegmann, K.; Olsen, L. F. Biochim. Biophys. Acta, in press. (4) Degn, H.; Olsen, L.; Perram, J. Ann. N.Y. Acad. Sci. 1979, 316, 623. (5) Olsen, L. Phys. Lett. 1983, 94A, 454. (6) Yokota, K.; Yamazaki, I. Biochemistry 1977, 16, 1913-1920. (7) Aguda, B.; Larter, R. J. Am. Chem. Soc. 1991, 113, 7913-7916. (8) Fed’kina, V. R.; Bronnikova, T. V.; Ataullakhanov, F. I. Biophysica 1992, 37, 781-789. (9) Olson, D.; Williksen, E.; Scheeline, A. J. Amer. Chem. Soc. 1995, 117, 2. (10) Bronnikova, T.; Fed’kina, V.; Schaffer, W.; Olsen, L. J. Phys. Chem. 1995, 99, 9309. (11) Hung, Y.-F.; Schreiber, I.; Ross, J. J. Phys. Chem. 1995, 99, 19801987. (12) Alexandre, S.; Dunford, B. Biophys. Chem. 1991, 40, 189. (13) Hauser, M. J. B.; Olsen, L. F. Biochim. Biophys. Acta, submitted.
JP9536110
