I. Bistability and Bursting Oscillations at Low Enzyme Concentrations

310

J. Phys. Chem. B 2001, 105, 310-321

Nonlinear Dynamics of the Peroxidase-Oxidase Reaction: I. Bistability and Bursting Oscillations at Low Enzyme Concentrations Tatiana V. Bronnikova* and William M. Schaffer Department of Ecology and EVolutionary Biology, The UniVersity of Arizona, Tucson, Arizona 85711

Lars F. Olsen Physical Biochemistry Group, Department of Biochemistry & Molecular Biology, SDU Odense UniVersity, Forskerparken 10, 5230, Odense M, Denmark ReceiVed: August 29, 2000; In Final Form: October 18, 2000

Under CSTR or semibatch conditions, the horseradish peroxidase (HRP)-catalyzed peroxidase-oxidase (PO) reaction evidences a wide range of nonlinear dynamical behaviors. Many of these regimes have proved to be predictable by a detailed model of the reaction first proposed in 1995. This model, which we refer to as BFSO, can also account for experimentally observed bifurcation sequences in response to varying concentrations of phenolic modifiers and rates of hydrogen donor input. Among those findings for which the model cannot account is the observation of bistability and bursting at low enzyme concentrations. This deficiency is important not only because these phenomena are biologically important but also because their existence requires a topology which, for the experimental circumstances in question, appears to be inconsistent with the model as originally formulated. In the present paper, we show that this deficiency can be remedied by the inclusion of an additional reaction whereby NADH and superoxide anion react in the presence of hydrogen ion to produce NAD radicals and hydrogen peroxide. Comparison of the modified model’s behavior with laboratory experiments suggests semiquantitative agreement between theory and observation. In particular, the model is able to reproduce experimentally observed responses to short-term perturbation by oxygen input suspension and the addition of hydrogen peroxide to the reaction mixture, as well as what was first described as “autonomous” switching between stable and oscillatory dynamics. Mathematically, addition of the new reaction makes possible the interaction of Hopf and hysteresis instabilities, as previously described in the BelousovZhabotinskyi reaction.

I. Introduction The first indications of chaos in a chemical oscillator were reported1 in 1977 for the peroxidase-oxidase (PO) reaction, a system in which oscillatory dynamics had previously been observed by Yamazaki and his associates2,3 in the mid-1960s. Since that time, the PO reaction has become a model system for the study of biochemical dynamics. Continuing interest in this reaction derives both from the extensive range of nonlinear behaviors evidenced in the laboratory and from the fact that peroxidase enzymes are widely distributed in the plant and animal kingdoms, where they function in a variety of critical capacities.4-7 There is also accumulating, albeit indirect, evidence that oscillations observed in vitro are relevant to the metabolism of living cells. In the case of horseradish peroxidase (HRP), the enzyme species most widely utilized experimentally, the traditional laboratory modifier, 2,4-dichlorophenol (DCP), can be replaced by natural compounds without compromising the system’s capacity for oscillatory dynamics.8,9 Other peroxidases have been shown to induce sustained oscillations as well,3,8,10,11 suggesting that the more extensive results for HRP may be broadly representative. The possibility of peroxidasemediated oscillations in vivo is further strengthened by the * Corresponding author. † The University of Arizona. ‡ SDU Odense University.

observation of damped oscillations in experiments utilizing bound enzyme12 and cell-free extracts obtained from horseradish roots.13 Although workers in the field often speak of the PO reaction, peroxidases actually catalyze the oxidation of a large number of hydrogen donors. In most cases, the oxidizing agent is H2O2, but for a few substrates, molecular oxygen can serve in its place. Such reactions are called peroxidase-oxidase reactions. In laboratory experiments, the hydrogen donor is often NADH, in which case the overall stoichiometry is

2NADH + O2 + 2H+ f 2NAD+ + 2H2O

(I)

Reaction I is typically studied under semibatch or CSTR conditions with continuing input of NADH and oxygen. Under these circumstances and in the presence of modifiers such as DCP and methylene blue (MB), the concentrations of O2, NADH, and various enzyme intermediates can oscillate.14-16 The reported oscillatory states include simple3 and mixedmode17,18,19 oscillations (MMOs), bursting oscillations,1,20 quasiperiodicity,17,21 and low-dimensional chaos22-24 of several types.25 In recent years, most of these behaviors have been shown to be reproducible by a mathematical model26 which implements the first 13 reactions in Table 1. This model, which we refer to

10.1021/jp003108+ CCC: $20.00 © 2001 American Chemical Society Published on Web 12/07/2000

Peroxidase-Oxidase Reaction

as BFSO, can also reproduce bifurcation sequences observed experimentally17,22,23,27 in response to variations in DCP concentration,26 the rate of NADH input,27,28 and the oxygen concentration in the gas phase.29,69 Not surprisingly, there are also experimental results which the model cannot explain. Of these, the observation of bistability and bursting oscillations1,20 at low enzyme concentrations is arguably the most important. Bistability and bursting are of widespread interest because of their demonstrable relevance to biological function.30-36 In addition, they also require a phase space topology not generated by the model for parameter values corresponding to the experimental conditions. It follows that for these conditions BFSO fails as a minimal model. This failure is not wholly unexpected. Schemes such as BFSO omit a large number of intermediate reactions (elementary steps) which, depending on the circumstances, can be expected to exert greater or lesser influence on reaction dynamics.6,16,37 The observation of model failure for low enzyme concentration raises two related questions:

J. Phys. Chem. B, Vol. 105, No. 1, 2001 311

(i) To what extent must the BFSO mechanism be modified, e.g., by incorporation of additional reactions, to make it compatible with experimental observations? (ii) To what extent, if any, do these modifications compromise previously published claims of model-data correspondence? In the present paper, we address the first question. Specifically, we observe that BFSO’s inability to induce bistability and bursting at low enzyme concentrations can be remediated by the inclusion of an additional reaction38,39 whereby NADH and superoxide anion, O2-, react in the presence of hydrogen ion to produce H2O2 and NAD radical, NAD•, i.e.,

NADH + H+ + O2- f NAD• + H2O2

(R14)

With regard to models of eq I, the history of this reaction is an on-off affair. Along with R15, R16, and several other elementary steps not included in Table 1, R14 was postulated by Yokota and Yamazaki40 and subsequently included in many of the socalled FAB3 models41-47 and their North American deriva-

312 J. Phys. Chem. B, Vol. 105, No. 1, 2001

Figure 1. Horseradish peroxidase (HRP) peroxidase-oxidase reaction, as modeled by the first 13 reactions in Table 1.

tives.48-50,70 More recently, Olson et al.16 argued against R14 and excluded it from their Urbanalator (U). This point of view was carried over to BFSO. Elsewhere,51 we observe that the addition of R14 does not invalidate previously published26,27,28 correspondences between BFSO and experiments carried out at [ET] g 0.9 µM. This answers the second question. II. Peroxidase-Oxidase Reaction Whereas the overall balance equations, e.g., eq I, of PO reactions are simple, the actual mechanisms (Figure 1) are rather complicated. Typically, they involve multiple enzyme states and the production of active species such as H2O2, NAD•, and superoxide anion, O2-. In the case of HRP, it is generally held7 that there are five principal enzyme species: compound I (coI or Per5+), compound II (coII or Per4+), ferric peroxidase (Per3+), ferrous peroxidase (Per2+), and compound III (coIII or Per6+), which differ according to oxidation state. The reactions by which these species are transformed into each other constitute a pair of oxidation-reduction loops, coupled by the successive reductions (R3 and R4) of coI and coII. Of utmost importance is the fact that the latter reactions generate NAD radicals via the oxidation of NADH. The first of the two feedback loops is the peroxidase cycle

Per3+ f coI f coII f Per3+ whereby Per3+ is oxidized (R2) to coI by H2O2. The second feedback loop, the oxidase cycle

Per3+ f coIII f coI f coII f Per3+ entails the production of coIII, which is then reduced to coI. The transformation of Per3+ to coIII can either be direct (R6), in which case the oxidant is O2-, or it can proceed via the formation of Per2+, i.e., Per3+ f Per2+ f coIII. In either case, a final reduction (R8) transforms coIII into coI. Additional complications to the oxidase cycle are discussed by Dunford.7 These include the monomolecular decomposition of coIII to coII (R15) or Per3+ (R16).7,52-55 Critical to the oxidase cycle are the NAD radicals which mediate the transformation of coIII to coI (R8) and, less importantly, that of Per3+ to Per2+ (R10). With each iteration of the cycle, two such radicals are produced by the reduction of coI to Per3+ (R3 and R4). Of these, one radical is consumed by the reduction of coIII to coI (R8) and a fraction of a radical by the reduction of Per3+ to Per2+.71 On the surface, these considerations might lead one to expect that NAD radicals would accumulate. In fact, not all the NAD radicals consumed

Bronnikova et al. participate in the reduction of enzyme intermediates. Some are converted to NAD+ when O2 transforms to O2- (R5) which, in turn, mediates the oxidation of Per3+ to coIII (R6). Still others combine (R9) to form NAD2 dimers which historically have been assumed6,16,37 to be unreactive. It follows that for a wide range of parameter values, the oxidase cycle is unable to sustain itself in isolation. Continuing inputs of NADH notwithstanding, the supply of NAD and superoxide anion radicals, along with that of H2O2, will often be exhausted, at which point the reaction shuts down. Herein lies the significance of the coupled cycles. Even when the oxidase cycle dominates (see below), it is the peroxidase cycle that sustains the reaction by utilizing H2O2 to replenish the supply of coI (R2) and thereby provide for the production of additional NAD radicals. While most of the reactions involving the various enzyme intermediates now seem to be well understood, many of the uncatalyzed reactions have not been resolved. In particular, reactions involving the modifiers, DCP and MB, need to be studied more carefully and their consequences to model dynamics assessed. In addition, the highly nonlinear character of the reaction network underlying (I) makes for the possibility that, at least under certain circumstances, seemingly insignificant reactions such as R14 can have important dynamical effects. In the present and in a forthcoming paper,51 we concentrate on this reaction which, as mentioned above, has been excluded from recent models because O2- is quite unreactive toward NADH.7,56 Actually, R14 consists of two reactions

HO2• + NADH f H2O2 + NAD•

(R14a)

and H+

O2- + NADH 98 H2O2 + NAD•

(R14b)

with k14a ) (1.8 ( 0.2 × 105) M-1 s-1 39 and k14b , 27 M-1 s-1.56 Since the pKa value for the dissociation of hydroperoxyl radical, HO2•, is 4.9, it follows that as the pH of the reaction mixture approaches this value, an increasing fraction of O2will exist in its protonated form. This accounts for the reported38,39 pH dependence of k14. Because of the rapid equilibrium between superoxide and the hydroperoxyl radical (HO2• a O2- + H+), we may consider reactions R14a and R14b as one reaction, R14, which proceeds at rate k14[O2-]T[NADH]. Here, [O2-]T ) [O2-] + [HO2•] is the total concentration of superoxide in solution, and

k14 )

Rk14a + k14b R+1

R ) 10(pKa-pH)

According to this, k14 should have a value of 7.0 × 104 M-1 s-1 at pH 5.1 and a value of 57 M-1 s-1 at pH 8.4. This does not square very well with the estimate of Land and Swallow of k14 , 27 M-1 s-1 at the latter pH. Hence, we suggest that the value reported by Nadezhdin and Dunford39 is an overestimate of k14a and, in the calculations that follow, assume a value approximately 2 orders of magnitude lower, which better matches the reported pH dependence of k14 observed by others.38,56 Another species which deserves more attention is NAD dimer. NAD2 was previously considered to be an inert species. However, recent experiments57 show that NAD2 reduces coI and coII with rate constants that are comparable to the rate constants for the reduction of these enzyme species by NADH. Furthermore, NAD2 was shown to reduce coIII to coI fairly rapidly, although the reported rate constant for this reaction57



Figure 2. Experimental evidence for bistability. Left: Transition from periodic behavior to a steady state induced by the addition (arrow) of H2O2. Right: Transition from a steady state to apparently periodic dynamics induced by the temporary shutdown (double arrow) of oxygen input. Reproduced with permission from Aguda et al.20

Figure 3. Spontaneous transitions from steady-state to oscillatory dynamics and back. Reproduced with permission from Aguda et al.20

is about two orders lower than the rate constant assumed for reaction R8 (see Table 1). It is also worth mentioning that reaction R8 has never been observed directly, and a radiolytic study at pH above 7.0 did not show a reaction between coIII and NAD•.58 We intend to address these problems in a forthcoming study. III. Experimental Observations of Bistability and Bursting Aguda et al.20 reported bistability (inducible switching from one apparently stable state to another) at low enzyme concentrations under acid conditions (pH 5.1) for various concentrations (20-40 µM) of DCP. This claim was based on experiments in which short-term perturbations induced switching between steady-state and periodic dynamics. The results of two such experiments are shown in Figure 2. At the left, a switch from periodic to steady-state behavior is induced by the addition of a small amount of H2O2. At the right, a transition from steadystate to periodic behavior is induced by briefly suspending the flow of oxygen into the reaction vessel. Figure 3 (also reproduced from their paper) illustrates what Aguda et al. described as “autonomous” switching between fixed-point and limit-cycle dynamics. As discussed below, these data may also be interpreted as multispike bursting oscillations, i.e., as resulting from motion on a single attractor. Bursting oscillations, but not bistability, have also been observed at low [ET] by Olsen and Degn, who reported 1-3 spike oscillations in their original paper.1 IV. Modeling the Reaction Efforts to model the PO reaction date back to Degn and Mayer59 who, together with Yokota and Yamazaki,40 Olsen and

Degn,14 and Fed’kina et al.,41 were the first to attempt to capture the dynamics of eq I mathematically. As discussed, for example, by Scheeline et al.,6 there are upwards of 35 possible intermediary reactions, some of which are listed in Table 1. This makes possible a large number of alternative mechanisms which differ with regard to the inclusion or omission of particular reactions. In the present paper, we focus on the model (BFSO) of Bronnikova et al.,26 which implements the first 13 reactions in Table 1. The result is a system of 10 differential equations with 16 parameters (Table 2). The latter include the 14 rate constants, k1 to k-13, the stock concentration of NADH, [NADH]st, and the equilibrium concentration of oxygen in the vessel, [O2]eq. The parameter values used in the present study are also given in Table 1. 1. BFSO Dynamics. Like U and in contrast to the FAB models,41-45,47 BFSO primarily engages the oxidase cycle. By this, we mean that coI is produced principally by the reduction of coIII rather than by the oxidation of Per3+. This is shown in Figure 4, wherein we plot the reaction rates, R2 ) k2[H2O2][Per3+] and R8 ) k8[coIII][NAD•], against time for k8 ) 1.37 × 108 M-1 s-1 and the remaining parameter values as in Table 1. Clearly, the area under R8 (coI produced by the oxidase cycle) greatly exceeds that under R2 (coI produced by the peroxidase cycle). Within the oxidase cycle itself, most coIII results from the direct oxidation (R6) of Per3+, with the consequence that only small amounts of Per2+ are produced and then principally at low values of k8. This latter observation may explain the conflicting experimental results regarding the presence of Per2+ during the reaction.60,61


Bronnikova et al. enzyme intermediates are driven to a steady state in which only Per3+ and coIII are present. This effect is reversible. Restoring k1 to its original value restarts the peroxidase cycle. This leads to coI formation and thence to the formation of NAD radicals. The latter allow for the resumption of coIII f coI conversion (R8) and the production of O2- (R5) which, in turn, permits the oxidation of Per3+ to coIII (R6). Not all models of the PO reaction include R1. In its absence, reaction shutdown can be avoided by the inclusion of reactions R14 (oxidation of NADH to NAD• by O2-), which occasions the production of H2O2 and R15 or R16 (decomposition of coIII to coII or Per3+), which respectively produces H2O2 or O2-.52-55 Like R1, these reactions provide for the replenishment of coI and the continuing production of NAD radicals. Indeed, summing reactions

NAD• + O2 f NAD+ + O2-

Figure 4. CoI production via reactions R2 and R8 in BFSO. k8)1.37 × 108 M-1 s-1.

(R5)

and

NADH + H+ + O2- f NAD• + H2O2

(R14)

yields

NADH + O2 + H+ f NAD+ + H2O2

(R1)

We emphasize that it is the substitution of R1 for R14 and R15 that distinguishes BFSO and its immediate antecedent, U, from models of the Russian school.41,43,45-47 The latter, in turn, derive from the original model of Yokota and Yamazaki,40 which included R1. The consequences of these variations in modeling strategy will be considered elsewhere. V. Simulations 1. Methods. The reactions in Table 1 were translated into differential equations, as indicated in Table 2. Terms shown in bold are those resulting from the addition of R14 to the original model. To document bistability, we modified eq II by observing that the concentrations of the various enzyme intermediates sum to a constant, [ET], and by replacing the differential equation, (d[Per3+]/dt) ) ‚‚‚, with the algebraic expression

[Per3+] ) [ET] - [coI] - [coII] - [coIII] - [Per2+] Figure 5. Shutdown of BFSO mechanism occasioned by the suppression of R1. The rate constant, k1, is set to zero at t ) 50 000 sec and restored to its original value at t ) 51 000 sec. During the shutdown, the concentrations of NADH and O2 rise, while those of the enzyme intermediates (not shown), CoI, CoII, and Per2+, drop to zero. Note the brevity of the shutdown relative to the time required for resumption of oscillations. [ET] ) 0.9 µM; k8 ) 1.37 × 108 M-1 s-1.

2. Importance of R1. Key to the BFSO mechanism is the nonenzymatic oxidation of NADH, reaction R1. This reaction

NADH + O2 + H+ f NAD+ + H2O2

(R1)

which was reintroduced by Olson et al.,16 generates hydrogen peroxide, thus driving the peroxidase cycle and thereby replenishing the supply of coI. Without R1, the BFSO mechanism comes to a halt for plausible parameter values. This is shown in Figure 5, wherein the effects of temporarily setting the rate constant, k1, to zero are illustrated. During the period of reaction shutdown, NADH accumulates, and [O2] rises to the maximum value possible. At the same time, the concentrations of all three active species, O2-, H2O2, and NAD•, drop to zero. As a result, both feedback loops are broken, and the concentrations of

(III)

This yields a “reduced” system consisting of nine ODEs. For the same parameter values and initial conditions, the two models, i.e., original and reduced, give equivalent results, as evidenced by test calculations. Numerical integration was accomplished using a doubleprecision implementation (LSODE) of Gear’s method62 for solving stiff systems of differential equations. Unless otherwise noted, parameter values are those given in Table 1. Of these, the rate constants, k1-k11, were obtained from the literature,5,22,23,26,39,51,54,56 while k12[NADH]stock, k13[O2]eq, and k-13 were chosen to match the experimental setup of Geest et al.22,23,72 Asymptotic trajectories for representative parameter values were computed by accumulating points at time intervals of 1 simulated second, following a transient of no less than 40 000 simulated seconds. Bifurcation diagrams were computed by recording the system’s asymptotic dynamics for 450+ values of a single parameter over a specified range. For each parameter value, a minimum of 150 successive maxima in simulated concentrations of oxygen were accumulated following a transient of 40 000100 000 simulated seconds. Except in cases where there was

Peroxidase-Oxidase Reaction reason to use the same initial conditions for different parameter values, the diagrams were computed in the usual way: For the first parameter value, a vector of starting concentrations was provided. For subsequent parameter values, the final concentrations for the preceding parameter value were used as initial conditions for the current parameter. Provided that the difference between successive parameter values is small, this procedure allows the computer to “shadow” the evolution of attractors obtained in response to a slowly varying parameter. CaVeat. In interpreting the bifurcation diagrams, the reader should bear in mind the following particulars: (i) Extracting a next amplitude (e.g., successive maxima) map from a time series is topologically equivalent to computing a Poincare´ section and consequently entails a reduction in dimension. Thus, simple cycles become fixed points; period-2 orbits, two-point cycles; tori, invariant loops, etc. (ii) This leads to certain ambiguities when map dynamics are projected onto a single dimension in diagrams which display a one-dimensional representation (the y-axis) of system dynamics for a range of parameter values. So, for example, regions of toroidal flow become indistinguishable from chaos; i.e., both correspond to dense bands in the diagram. (iii) An additional complication regards the representation of fixed-point behavior. In this case, there being no extrema, the only sensible course is to plot the fixed points themselves. As a result, a curve on the bifurcation diagram can correspond either to a sequence of period-1 cycles, in which case the cycles’ maxima are viewed, or to a sequence of equilibria, in which case the actual fixed points are viewed. (iv) In light of the cases studied below, it is instructive to consider the form of the bifurcation diagram for the case of a cycle which emerges smoothly (i.e., via a supercritical Hopf bifurcation) from a stable point. A three-dimensional plot of the continuous dynamics would show a paraboloid (the cycles) emerging from a point. Conversely, a two-dimensional diagram, in which cyclic maxima are plotted against the bifurcation parameter, would evidence only a change in slope of a continuous curve. In what follows, we shall call attention to such matters when the nature of the dynamics cannot be determined from a onedimensional projection. 2. Dynamical Consequences of Adding R14 to the Model. As parametrized in Table 1, BFSO is incapable of generating either bistability or bursting. Moreover, in all instances investigated, reducing total the enzyme concentration leads to nothing more interesting than low-amplitude cycles of period-1. We illustrate this observation in Figure 6 (left). Here, we display a series of bifurcation diagrams in which the total enzyme concentration, [ET], is varied from 0.10 to 1.2 µM. The several diagrams correspond to decreasing values (top to bottom) of the rate constant, k8, which was used previously26 as a proxy variable for [DCP].73 In each case, decreasing [ET] from 0.9 µM (the value used to simulate the DCP experiments) leads to an overall reduction of dynamical complexity, if such exists. In all cases, the dynamics at low [ET] are simple cycles. Periodic Dynamics and Bursting at Low [ET]. Recalculating these diagrams with R14 added to the reaction mechanism and k14 ) 3 × 102 M-1 s-1 yields a very different picture, as shown in Figure 6 (right). In particular, we note the following: (i) For k8 ) 2.55 × 108 M-1 s-1, there is an abrupt change from period-1 dynamics to a steady state at low (∼0.25 µM) [ET]. As discussed above, both behaviors manifest themselves in bifurcation diagrams as one-dimensional curves. In the present case, the curve occupying most of the diagram, i.e., the one

J. Phys. Chem. B, Vol. 105, No. 1, 2001 315 with the “bubble,” corresponds to successive maxima of a period-1 (period-2 in the bubble) cycle, while the much shorter curve below and to the left is a curve of equilibria. (ii) For k8 on the interval, [1.05 × 108, 2.25 × 108 M-1 s-1], a region of complex behavior exists at low [ET]. (iii) For k8 ) 1.05 × 108 M-1 s-1 and k8 ) 2.25 × 108 M-1 s-1, the region of complex dynamics terminates (at the left) before the switch to steady-state behavior. For the remaining cases, i.e., k8 ) 1.35 × 108 M-1 s-1, 1.65 × 108 M-1 s-1, and 1.95 ×108 M-1 s-1, complex dynamics persist up to the transition to equilibrium. The complex dynamics at low [ET] turn out to be quasiperiodic. This is shown in Figure 7, wherein we display representative time series and next-amplitude maps. The former evidence modulated periodicity, i.e., the addition of a second frequency, while the latter reveal the characteristic invariant loops of tori in section. With decreasing [ET] (and fixed k8), quasiperiodicity arises at the right (high [ET]) by torus (secondary Hopf) bifurcation and is lost at the left (low [ET]) either by reverse torus bifurcation (k8 ) 1.05 × 108 and 2.25 × 108 M-1 s-1) or interior crisis (k8 ) [1.35×-1.95 × 108 M-1 s-1). In the former case, further [ET] reduction occasions a transition from period-1 oscillations to steady-state behavior. In the case of torus death by interior crisis, the corresponding transition is from quasiperiodicity to a steady state. At low [ET], the tori assume a distinctly heteroclinic character (Figure 7, bottom). In the time domain, this behavior manifests itself as bursting oscillations which bear a striking resemblance to the “autonomous,” i.e., self-induced, switching (Figure 3) reported by Aguda et al.20 Bistability. The diagrams in Figure 6 were computed using eq II, (i.e., 10 state variables) with the initial concentration of Per3+ on the interval [0.1, 1.2 µM] and with [O2] ) 16.7 µM and the remaining initial concentrations set to zero. To visualize bistability, we use the reduced eq III with [ET] as the bifurcation parameter. Computing forward and backward diagrams reveals the existence of a region of bistability in which a stable fixed point coexists with a stable cycle or torus. A representative example is shown in Figure 8. Mathematical Basis. The aforementioned transitions appear to be the consequence of a codimension-3 bifurcation involving the interaction of Hopf and hysteresis instabilities,63,64 as discussed, for example, by Richetti et al.65 and Barkley et al.66 in the case of the Belousov-Zhabotinskii reaction. A more mundane representation of the dynamical essentials is attempted in Figure 9. Here we display pairs of three-dimensional projections of the phase space. Each pair of projections illustrates the evolution of trajectories based on two initial conditions: [NADH](0) ) 150.0 µM (left) and [NADH](0) ) 0 µM (right) for a single value of [ET] with the remaining initial concentrations equal to zero. As one moves down the figure, the total enzyme concentration is decreased from 0.45 to 0.25 µM. Inspection of the paired trajectories suggests that there are three equilibria, all of which are initially (i.e., at [ET] ) 0.45 µM) saddle foci. Of these, the left- and right-most equilibria have one-dimensional outsets (unstable manifolds), while the middle saddle has a two-dimensional outset. In addition, we note the following: (i) Figure 9a. At [ET] ) 0.45 µM, the trajectories based on both initial conditions are attracted to a periodic orbit close to the high NADH saddle. (ii) Figure 9b. As [ET] is reduced, this orbit turns into a torus to which trajectories based on both initial conditions are


Bronnikova et al.

Figure 6. Effects of varying initial concentration of Per3+ in the original BFSO model (left) and with the addition of R14 (right). For each parameter value, the starting concentrations, [CoI], [CoII], [CoIII], and [Per2+], of the remaining enzyme species all equal zero. From top to bottom: (a) k8 ) 2.55, 2.25, and 1.95 × 108 M-1 s-1 and (b) k8 ) 1.65, 1.35, and 1.05 × 108 M-1 s-1.


Figure 7. Periodic (top) and quasiperiodic (middle and bottom) behavior induced by BFSO with the addition of R14 at low [ET]. Left: [O2] time series. Right: Next amplitude maps constructed by plotting successive maxima in [O2]. Top to bottom: [ET] ) 0.45, 0.36, and 0.32 µM. k8 ) 2.0 × 108 M-1 s-1.

Figure 8. Bistability observed in response to decreasing (left) and increasing (right) total enzyme concentration, [ET], in BFSO supplemented by R14. k8 ) 2.20 × 108 M-1 s-1; k14 ) 3.0 × 102 M-1 s-1.

attracted. (iii) Figure 9c. With further reductions in [ET], the torus grows and acquires a distinctly heteroclinic character. At this point, it is still the sole attractor. (iv) Figure 9d. With additional reductions in [ET], the middle fixed point is stabilized. At this point, heteroclinic orbits connecting it and remaining saddles come into existence. Now there are two attractors. Trajectories based on [NADH](0) ) 0 tend to the newly stabilized fixed point, while trajectories based on [NADH](0) ) 150 µM still tend to the torus. (v) Figure 9e. Still further reductions in [ET] lead to destruction of the torus, which dies when it collides with the point attractor’s basin of attraction. With destruction of the torus, trajectories based on both initial conditions tend to a steady state, and all oscillatory tendencies are lost. The trajectories in Figure 9 were computed for k8 ) 2.2 × 108 M-1 s-1, in which case torus destruction at low [ET] results from interior crisis. For other values of k8, torus destruction results from a secondary Hopf bifurcation, in which case it is a

J. Phys. Chem. B, Vol. 105, No. 1, 2001 317 periodic orbit, as opposed to a torus, which collides with the basin boundary of the fixed point. Dependence on k14. A deeper understanding of the dynamics is obtained if the loci of changes in behavior in the k14-[ET] plane are plotted. The result is a pair of curves, as shown in Figure 10. One of these (O-O-O) is a curve of Hopf bifurcations. For parameter pairs corresponding to points in the plane above this curve, the system evidences fixed-point behavior. This corresponds to those cases in Figure 9 in which only the middle fixed point is an attractor. For parameter values below this curve, periodic (or quasiperiodic) dynamics are manifest. Emerging from the curve of Hopf bifurcations is a curve of cyclic folds (+ + +) where a pair of periodic orbits (one stable, the other not) come into existence. To the right of this curve, periodic (or quasiperiodic) behavior exists in the phase space. Between the two curves, this behavior coexists with fixed-point dynamics. This is the region of bistability. Critical to the existence of bistability is the fact that what happens when one crosses the curve of Hopf bifurcations depends on whether one is to the right or to the left of the fold curve. If one crosses to the right, the unstable periodic orbit created by the fold bifurcation is absorbed by the middle fixed point, which is thereby destabilized. In this case, the Hopf bifurcation is subcritical and can induce an abrupt change in system dynamics, i.e., from fixed point to oscillatory dynamics. If, on the other hand, one crosses the Hopf curve to the left of the fold curve, the Hopf bifurcation will be supercritical and a stable periodic orbit and will emerge smoothly from the middle fixed point. In short, bistability requires experimental circumstances corresponding to points in the k14-[ET] plane, which lie between the two curves. For the parameter region studied here, the curve of Hopf bifurcations is bounded away from k14 ) 0. Consequently, the codimension-2 bifurcation point67 at which the curve of fold bifurcations emerges from the Hopf curve is also bounded away from k14 = 0. This is the basis for our ascertion that BFSO as originally formulated is inconsistent with experimental observations of bistability and bursting at low [ET]. 3. Simulating Experimental Observations of Bistability. Having demonstrated the existence of bistability and bursting in BFSO supplemented by R14, we next inquire as to whether transitions from steady-state to oscillatory (periodic or quasiperiodic) dynamics and back can be induced in simulation by temporary O2 shutdown and the addition of H2O2. Figure 11 answers this question in the affirmative. Here, we show a single trajectory (top) subject to “spiking” with H2O2 and temporary O2 input shutdown. As observed experimentally, the former induces a transition from periodic to steady-state dynamics, whereas the latter perturbation induces the reverse transition. Magnifications (bottom, left, and right) of the transitions bear a remarkable resemblance to the experimental results20 reproduced in Figure 2. VI. Discussion How well does the theory square with observation? On the plus side, the R14 extended version of BFSO (dare we call it BFSO-14) reproduces the following experimental findings: (i) Bistability at low [ET] for a wide range of DCP concentrations. (ii) Inducible transitions between steady-state and oscillatory dynamics in response to short-term oxygen input shutdown and the addition of H2O2. (iii) “Autonomous”, i.e., self-induced, transitions from steadystate to oscillatory dynamics.


Bronnikova et al.



Figure 9. Effect of decreasing enzyme concentration on the time evolution of simulated systems based on different initial conditions: [NADH](0) ) 150 (left) and 0 (right) µM, with [O2](0) ) 16.7 µM and the initial concentrations of the remaining reactants set to 0. k8 ) 2.0 × 108 M-1 s-1; k14 ) 3 × 102 M-1 s-1. Total enzyme concentrations as follows: (a) 0.45 µM; (b) 0.41 µM; (c) 0.365 µM; (d) 0.275 µM; (e) 0.25 µM.

Figure 10. Dynamics in the k8-[ET] parameter plane for BFSO with the addition of R14. Emerging from a curve of Hopf bifurcations (OO-O) is a curve (+ + +) of cyclic folds. The former corresponds to a transition from steady-state to fixed point behavior, while the latter delimits the region of bistability. See text for additional details.

With regard to bistability, the model predicts coexisting attractors at low [ET] (as observed) for a wide range of k8 values, which we interpret to mean a wide range of DCP concentrations (also observed). Of the two attractors, one is always a steady state; the other is a periodic orbit or quasiperiodic motion on a torus. As a practical matter, the latter will often, i.e., in the absence of pronounced heteroclinicity, be experimentally indistinguishable from noise-corrupted periodicity. Significant, we believe, is the fact that the predicted steady-state value of [O2] lies between the high and low values of [O2] evidenced by the coexisting periodic (quasiperiodic) attractor (Figures 9 and 11). This is in agreement with the experimental observations and in disagreement with the predictions of an earlier model.48 In addition, the cycle amplitude and period are roughly comparable to the experimental observations (compare Figures 2 and 11). Finally, the duration of the O2 input shutdown (90 s) required to effect a transition of steady-state to oscillatory behavior in simulation is within the range (30-90 s) employed experimentally.

Figure 11. Simulated (BFSO with addition of R14) transitions from periodic to steady-state and back, induced respectively by the addition of 20 µM H2O2 (t ) 30,000 s) and the suspension of O2 input for 90 simulated seconds (t ) 60,000 s). Top: The full time series. Bottom: Magnification of the transitions. Left: Periodic to steady state transition induced by H2O2 addition. Right: Steady state to periodic transition. Compare with Figure 2.

There are also discrepancies between model predictions and the observations. A minor problem is that the enzyme concentrations (0.25-0.35 µM) for which the model predicts bistability is somewhat lower than that (0.4-0.5 µM) observed experimentally. However, this is specific to setting k14 ) 3 × 102 M-1 s-1. For higher values of the rate constant, the values of [ET] at which bistability and bursting are observed increase. Of greater concern is the amount of H2O2 required to quench oscillatory dynamics, which, in simulation, is 3-4 time larger than the quantity used experimentally. Mitigating this divergence is the fact that addition of H2O2 induced dynamical change in only about half of the experiments. In addition, smaller amounts of H2O2 can suffice in simulation, but in this case, the transition is more gradual than what is observed in the lab. Nonetheless, this discrepancy reminds us that BFSO-14 is not the final word in minimal models of reaction I.

320 J. Phys. Chem. B, Vol. 105, No. 1, 2001 As emphasized by Scheeline and his associates,6,16,37 models such as BFSO implement only a fraction of the possible elementary steps in reaction I. The omissions include the monomolecular decomposition of coIII (R15, R16), reactions involving NAD2,57 and explicit hypotheses regarding the action of MB and DCP. Rather than attempting to deal with the enormous complexity that would result from the inclusion of all of these reactions, the present paper pursues a more conservative approach, which is to determine the ways in which a promising model must be modified to allow it to account for data which is at variance with its predictions. In short, we have elected to focus on a minimal model, as opposed to a complete one. Elsewhere,51 we point out that what constitutes a minimal model depends on the data which one seeks to explain. In this spirit, we would argue that the Urbanalator (U), from which BFSO derives, is a satisfactory minimal model for experiments that place the system in regions of the parameter space of a hypothetical “complete” model that correspond to period-1 dynamics. For experiments, placing the system in regions of the full space corresponding to complex dynamics, U, by virtue of its apparent inability to generate anything other than simple cycles, must be replaced by BFSO or some other scheme, e.g., one of the FAB models, compatible with greater dynamical complexity. Correspondingly, in the case of experiments evidencing bistability and bursting at low enzyme concentrations, BFSO must itself be replaced by yet another model, which, according to the present suggestion, is BFSO with the addition of R14. Of course, phrases such as “minimal model” and “which can account for the experimental observations” are loaded. One possible criterion, certainly not the only one, is that a model should be able to generate dynamics which are topologically equivalent to what is observed experimentally. In fact, topological equivalence is too strict a requirement. What one really wants is equivalence with regard to features that are experimentally resolvable.51 By this standard, models such as U and BFSO, in which the regeneration of coI requires the production of H2O2 via by R1 lack the ability to generate the requisite topology for bistability and bursting at low [ET]. Indeed, most of the published BFSO bifurcation sequences are nothing more than complicated examples of period-bubbling. It is therefore both interesting and relevant to experimental observation that the addition of R14 produces a new topology at low [ET]. The decision to proceed in a stepwise manner was further motivated by a desire to understand mathematically how variations in chemistry translate into variations in dynamics. Such an approach allows more complex schemes to be formulated as perturbations of simpler ones. This is a timehonored tradition in mathematical modeling which, on occasion, can lead68 to deep insights regarding the origins of experimentally observed phenomenology. Indeed, Figure 10 summarizes the results of just such an analysis, i.e., with k14 ) 0, we are back to BFSO in its first incarnation. The fact, as already commented upon, that the topology underlying bistability and bursting does not arise until k14 > 0 is key to all the results reported in this communication. In conclusion, it worth pointing out that our observations are very much of a piece with those of Fed’kina and Bronnikova, who observed that the addition of periodic forcing or R15 to a FAB model lacking R1 induces a topological change similar to that described here. In short, the comparative analysis of alternative models of reaction I as pioneered by the Russian school, but more closely tied to experimental observation, would

Bronnikova et al. appear to be an important next step in the quest for a definitive model of PO dynamics in the laboratory. Acknowledgment. L.F.O. is supported by a grant from the Danish Natural Science Research Council. Note Added after ASAP Posting This article was posted ASAP with an error in an equation on 12/07/00. The corrected version was posted on 12/15/00. References and Notes (1) Olsen, L. F.; Degn, H. Nature 1977, 267, 177-178. (2) Yamazaki, I.; Yokota, K.; Nakajima, R. Biochem. Biophys. Res. Commun. 1965, 21, 582-586. (3) Nakamura, S.; Yokota, K.; Yamazaki, I. Nature 1969, 222, 794. (4) Epstein, I. R. Phys. D 1991 51, 152-160. (5) Everse, J.; Everse, K. E.; Grisham, M. B. Peroxidases in Chemistry and Biology; CRC Press: Boca Raton, FL, 1991; Vols. I and II. (6) Scheeline, A.; Olson, D. E.; Williksen, E. P.; Horras, G.; Klein, M. L.; Larter, R. Chem. ReV. 1997, 97, 739-756. (7) Dunford, H. B. Heme Peroxidases; John Wiley: New York, 1999. (8) Kummer, U.; Hauser, M. J. B.; Wegmann, K.; Olsen, L. F.; Baier, G.. J. Am. Chem. Soc. 1997, 119, 2084-2087. (9) Hauser, M. J. B.; Olsen, L. F. Biochemistry 1998, 37, 2458-2469. (10) Kummer, U.; Valeur, K. R.; Baier, G.; Wegmann, K.; Olsen, L. F. Biochim. Biophys. Acta 1996, 1289, 397-403. (11) Valeur, K. R.; Olsen, L. F. Biochim. Biophys. Acta 1996, 1289, 377-384. (12) Cook, L.; Larter, R.; Shen, P.; Geest, T. J. Phys. Chem. 1993, 97, 9060-9063. (13) Møller, A. C.; Hauser, M. J. B.; Olsen, L. F. Biophys. Chem. 1998 72, 63-72. (14) Olsen, L. F.; Degn, H. Biochim. Biophys. Acta 1978, 523, 321334. (15) Olsen, L. F. Phys. Lett. A 1983, 94, 454-457. (16) Olson, D. L.; Williksen, E. P.; Scheeline, A. J. Am. Chem. Soc. 1995, 117, 2-15. (17) Hauck, T.; Schneider, F. W. J. Chem. Phys. 1993, 97, 391-397. (18) Hauck, T.; Schneider, F. W. J. Chem. Phys. 1994, 98, 2072-2080. (19) Hauser, M. J. B.; Andersen, S.; Olsen, L. F. In IV International Symposium of Plant Peroxidases; Obinger, C., Burner, U., Ebermann, R., Penel, C., Greppin, H., Eds.; Rochat-Baumann: Gene`ve, Switzerland, 1996; pp 88-93. (20) Aguda, B. D.; Frisch, L. L. H.; Olsen, L. F. J. Am. Chem. Soc. 1990, 112, 6652-6656 1990. (21) Samples, M. S.; Hung, Y.-F.; Ross, J. J. Phys. Chem. 1992, 96, 7338-7343. (22) Geest, T.; Steinmetz, C. G.; Larter, R.; Olsen, L. F. J. Phys. Chem. 1992, 96, 5678-5680. (23) Steinmetz, C. G.; Geest, T.; Larter, R. J. Phys. Chem. 1993, 97, 5649-5653. (24) Hauser, M. J. B.; Olsen, L. F. J. Chem. Soc., Faraday Trans. 1996, 92, 2857-2863. (25) Ro¨ssler, O. E. Ann. N.Y. Acad. Sci. 1979, 316, 376-392. (26) Bronnikova, T. V.; Fed′kina, V. R.; Schaffer, W. M.; Olsen, L. F. J. Phys. Chem. 1995, 99, 9309-9312. (27) Hauser, M. J. B.; Olsen, L. F.; Bronnikova, T. V.; Schaffer, W. M. J. Phys. Chem. B 1997 101, 5075-5083. (28) Bronnikova, T. V.; Schaffer, W. M.; Hauser, M. J. B.; Olsen, L. F. J. Phys. Chem. B 1998 102, 632-640. (29) Bronnikova, T. V.; Schaffer, W. M.; Olsen, L. F. J. Chem. Phys. 1996, 105, 10849-10859. (30) Dean, P. M.; Matthews, E. K. J. Physiol. (London) 1970, 210, 255264. (31) Hudson, J. L.; Hart, M.; Marinko, D. J. Chem. Phys. 1979 71, 1601-1606. (32) Deschenes, M.; Roy, J. P.; Steriade, M. Brain Res. 1982, 239, 289293. (33) Decroly, O.; Goldbeter, A. J. Theor. Biol. 1987, 124, 219-250. (34) Chay, T. R.; Cook, D. L. Math. Biosci. 1988, 90, 139-153. (35) Ashcroft, F.; Rorsman, P. Prog. Biophys. Mol. Biol. 1989, 54, 87143. (36) Holden, L.; Erneux, T. J. Math. Biol. 1993, 31, 351-365. (37) Kirkor, E. S., Scheeline, A.; Hauser, M. J. B. Anal. Chem. 2000, 72, 1381-1388. (38) Olsen, L. F. Biochim. Biophys. Acta 1978, 527, 212-220. (39) Nadezhdin, A.; Dunford, H. B. J. Phys. Chem. 1979, 83, 19571961. (40) Yokota, K.; Yamazaki, I. Biochemistry 1977, 16, 1913-1918.

Peroxidase-Oxidase Reaction (41) Fed’kina, V. R.; Ataullakhanov, F. I.; Bronnikova, T. V.; Balabaev, N. K. Stud. Biophys. 1978, 72, 195-202. (42) Fed’kina, V. R.; Ataullakhanov, F. I.; Bronnikova, T. V. Biophys. Chem. 1984, 19, 259-264. (43) Fed’kina, V. R.; Bronnikova, T. V. In Sinergetics-86; Shtinnitza: Kishinev, USSR, 1986; pp 126-127. (44) Fed’kina, V. R.; Ataullakhanov, F. I.; Bronnikova, T. V. Theor. Exp. Chim. 1988, 24, 172-178. (45) Fed’kina, V. R.; Bronnikova; V. R.; Ataullakhanov F. I. Biofizika 1992, 37, 885-892. (46) Fed’kina, V. R.; Bronnikova, T. V. Biofizika 1994, 39, 599-607. (47) Fed’kina, V. R.; Bronnikova, T. V. Biofizika 1995, 40, 36-47 1995. (48) Aguda, B. D.; Larter, R. J. Am. Chem. Soc. 1990, 112, 21672174. (49) Aguda, B. D.; Larter, R. J. Am. Chem. Soc. 1991, 113, 79137916. (50) Hung, Yu-F.; Schreiber, I.; Ross, J. J. Phys. Chem. 1995, 99, 19801987. (51) Schaffer, W. M.; Bronnikova, T. V.; Olsen, L. F. To be submitted. (52) Yamazaki, I.; Yokota, K.; Tamura, M. In; The Chemistry of Hemes and Hemoproteins; Chance, B., Estabrook, R. W., Yonetani, T., Eds.; AP: New York, 1966. (53) Tamura, M.; Yamazaki, I. Biochem. (Tokyo) 1972, 71, 311-319. (54) Metodiewa, D.; Dunford, H. B. Arch. Biochem. Biophys. 1989, 272, 245-250. (55) Adediran, S. A. Arch. Biochem. Biophys. 1996, 327, 279-285. (56) Land, E. J.; Swallow, A. J. Biochim. Biophys. Acta 1971, 234, 3442. (57) Kirkor, E. S.; Scheeline, A. Eur. J. Biochem. 2000, 267, 50145022. (58) Kobayashi, K.; Hayashi, K.; Swallow, A. J. Biochemistry 1990, 29, 2080-2084 1990. (59) Degn, H.; Mayer, D. Biochim. Biophys. Acta 1969, 180, 291-301. (60) Hauser, M. J. B.; Olsen, L. F., In IV International Symposium of Plant Peroxidases; Obinger, C., Burner, U., Ebermann, R., Penel, C., Greppin, H., Eds.); Rochat-Baumann: Gene`ve, Switzerland, 1996; pp 8287.

J. Phys. Chem. B, Vol. 105, No. 1, 2001 321 (61) Hauser, M. J. B.; Lunding, A.; Olsen, L. F. Phys. Chem. Chem. Phys. 2000, 2, 1685-1692. (62) Hindmarsh, A. C. ACM Signum Newslett. 1980, 15, 10. (63) Guckenheimer J.; Holmes, P. Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields; Springer-Verlag: New York, 1983. (64) Kuznetsov, Y. A. Elements of Applied Bifurcation Theory; SpringerVerlag: New York, 1995. (65) Richetti, P.; Roux, J. C.; Argoul, F.; Arneodo, A. J. Chem. Phys. 1987, 86, 3339-3356. (66) Barkley, D.; Ringland, J.; Turner, J. S. J. Chem. Phys. 1987, 87, 3812-3820. (67) Khibnik, A. I.; Kuznetsov, Yu. A.; Levitin, V. V.; Nikolaev, E. V. InteractiVe LOCal BIFurcation Analyzer, personal communication, 19901992. (68) King, A. A.; Schaffer, W. M. J. Math. Biol. 1999, 39, 439-469. (69) The name of BFSO is an acronym reflecting the workers (Bronnikova, Fed’kina, Schaffer, and Olsen) who proposed the model in 1995.26 (70) Like BFSO, the appellation, FAB, is an acronym reflecting the names (Fed’kina, Ataullakhanov, and Bronnikova) of those workers responsible for a series of models published principally in the Russian Journal of Biophysics in the years 1978-1995. (71) Fractional consumption results from the fact that some, generally most, of the oxidation of Per3+ to coIII proceeds via R6, in which case the sole oxidant is O2-. (72) The values for k1-k12 reported in the literature and the deviations therefrom are given by Bronnikova et al.;26 for k14, by Nadezhdin and Dunford39 and Land and Swallow51 (see below also); for k15, by Tamura and Yamazaki;54 and for k16, by Metodiewa and Dunford.55 The constants, k13 and k-13 are, of course, specific to the experimental setup. (73) Use of proxy variables for [DCP] and [MB] is necessitated by the fact that BFSO does not include the concentrations of these modifiers as explicit variables.

I. Bistability and Bursting Oscillations at Low Enzyme Concentrations

Recommend Documents