The peroxidase-oxidase (PO) reaction is an important example of how oscillating reactions arise in living organisms. The reaction is important mathematically because it exhibits many characteristics of chaos. Of these characteristics, sensitive dependence on initial conditions illustrates how predicting the future state of the PO reaction is nearly impossible. For this reason, chaotic behavior in living organisms presents many obstacles to chemists and biologists in their attempts to predict how systems will react to perturbations. This paper explores the chaotic behaviors of the PO reaction, using a system of four differential equations as a model. The topics analyzed are timeseries data, chaotic attractors, bifurcations, tests for sensitive dependence, Lyapunov exponents, and one-dimensional time delay embedding.
What is Chaos?
Many everyday occurrences act in unpredictable ways. Weather, the stock market, billiards, and planetary orbits are all examples of systems that have a degree of uncertainty due to either a myriad of factors influencing them or some other complexity. Chaos and the study of dynamical systems provide tools to model and analyze these complex and often seemingly random systems. Sensitive dependence on initial conditions is a classical characteristic of chaos. A system has sensitive dependence if a small perturbation in the initial conditions of the system results in rapidly increasing deviations in the state of the system. A common example is the butterfly effect on weather: the small perturbations in the wind from a butterfly flapping its wings in Asia could build up to cause a hurricane in the southeastern United States. A good way to test for sensitive dependence is using Lyapunov exponents. Lyapunov exponents measure the change in the natural log of the difference between the perturbed and the unperturbed state of a system at a given time. Thus a Lyapunov exponent greater than zero is characteristic of chaos.
Another characteristic often observed in chaotic systems is a period doubling route to chaos. For timeseries data, as in this paper, period doubling refers to a doubling of the number of maxima (as one of the parameters is varied monotonically) in the timeseries data. The number of maxima continues to double until the chaotic region is reached, in which there are infinite maxima. Period doubling is most easily illustrated in bifurcation diagrams. The bifurcation diagrams show the maxima of the timeseries data on the vertical axis and the value of the parameter being varied on the horizontal axis. The location that the lines split on these diagrams indicate where period doubling occurs. For a good introductory text on chaos, including more formal definitions of the terms above, see (1).
Oscillating Reactions
Less than a century ago, oscillating chemical reactions were thought to be nothing more than erroneous results caused by chemical impurities. Few believed that a reaction, in proceeding to equilibrium, could have intermediates with oscillating magnitudes of concentration. The Lotka-Volterra Oscillator and the Belousov-Zhabotinsky (BZ) reactions were two of the first oscillating reactions to be studied extensively and universally accepted to have oscillating concentrations of reactants (1). Finding chaos in oscillating reactions is a much more recent undertaking, gaining popularity with the rise of computers. In, fact there is still debate as to whether oscillating chemical reactions truly exhibit chaos, or if they are merely oscillatory with very large period (2). However since the discovery and eventual acceptance of oscillatory reactions, many more such reactions have been designed and discovered. One of the main areas of discovery of oscillating reactions is in living systems. The peroxidase-oxidase (PO) reaction is a prime example of such a reaction. It was one of the first reactions outside of the BZ reactions to be classified as oscillating and chaotic, and is an extensively studied example of an in vivo oscillating reaction (the reaction can be carried out both in vivo and in vitro). Molecular oscillating reactions are an important area of study, as they are essential to understanding the more complex oscillating systems of organism (e.g. a heartbeat) (1).
The Peroxidase-Oxidase Reaction
The peroxidase-oxidase reaction is an enzyme-catalyzed redox reaction. Nicotinamide adenine dinucleotide (NADH) is oxidized and molecular oxygen acts as an electron receiver. The net reaction is
2NADH + O2 + 2H+ g 2NAD+ + 2H2O
O2 and NADH are continuously replenished and products are continuously removed from the experimental system via the use of a continuous-flow, stirred tank reactor (CSTR). Experimentation has shown that the PO reaction exhibits oscillatory behavior and a period-doubling route to chaos (3). These characteristics have been effectively modeled using a simplified eight-step mechanism. The Olsen model (4) is
k1 k2 k3
B + X g 2X 2X g 2Y A + B + Y g 3X
k4 k5 k6
X g P Y g Q X0 g X
k7 k8
A0 g A B0 g B
where A and B are reactants (O2 and NADH respectively), P and Q are the products, and X and Y are reaction intermediates. Experimentally, [X] corresponds with [NAD] and [Y] corresponds with the concentration of oxyferrous peroxidase (compound III), however it is still unknown how well these intermediate variables correlate with the true intermediates. It is also important to note that k1 = [enzyme], the concentration of peroxidase enzyme, and there is also a strong correlation between k3 and the concentration of 2,4-dichlorophenol, [DCP] (3). Thus k1 and k3 are variable parameters. Oscillatory behavior arises because of competing mechanisms in the PO reaction, appearing in the Olsen model as follows: call mechanism I the net reaction of the autocatalytic production of X from B and X, and the production of 2Y from 2X.
Mechanism I
B + X g 2X
2X g 2Y
B + X g 2Y
Mechanism I dominates the PO reaction when the concentration of X is high, but uses up X to create Y. The rate of mechanism I can be varied by changing the concentration of peroxidase enzyme, k1. Once the concentration of X falls below some critical value,
[X]crit, a second mechanism, mechanism II, takes over. Mechanism II is the termolecular reaction of A, B, and Y to form X.
Mechanism II
A + B + Y g 3X
Mechanism II dominates when the concentration of Y is high, using up Y and turning it to X. Once the concentration of Y falls below some [Y]crit, mechanism I takes over once again. Thus the concentrations of the intermediates X and Y oscillate. Since the rate at which the reactants are converted into intermediates depends on the concentration of the intermediates, the concentrations of the reactants also oscillate and are strongly correlated. See figure 1 for an example of the correlation of [O2] and [NADH]. Also note that in mechanism II, the reaction rate is governed by k3, the concentration of DCP. Since k1 and k3 are simple to control experimentally, there is an easy way to change the rate constants that govern mechanisms I and II. In fact k1 and k3 determine the characteristics of how the concentrations of A, B, X, and Y vary with time, and if chaos arises. In fact, either decreasing k1 or increasing k3 results in a period-doubling route to chaos (3).
Mathematically Approximating the Olsen Model
Methods
All numerical analyses were carried out using MATLAB version 7.2.0.232 (R2006a). Systems of differential equations were solved using ode45 (a Runge-Kutta numerical approximation method). Timeseries data were numerically compared by fitting a cubic spline (using the MATLAB ‘fit’ toolbox) to the data from ode45 and then interpolating values from this fit. Linear regressions were also carried out using the ’fit’ toolbox. The Olsen model can be approximated by the following system of four first-order differential equations:
A = k7(A0-A) – k3ABY
B = k8 – k1BX – k3ABY
X = k1BX – 2k2X2 + 3k3ABY – k4X + k6
Y = 2k2X2 – k5Y – k3ABY
The following parameters were used in all of the calculations in this paper: k2 = 250, k4 = 20, k5 = 5.35, k6 = 10-5, k7 = 0.1, k8 = 0.825, and A0 = 8. Values of k1 and k3 are indicated in the diagrams. Initial conditions, except where indicated, are [O2] = A = 6, [NADH] = B = 58, [NAD] = X = 0, and [Co III] = Y = 0. Phase space diagrams were computed to compare [O2], [NADH], [NAD.], and [Co III]. No matter what values of k1 and k3 were used to compute the phase space, they all resulted in attractors. Chaotic attractors (as in figure 2) resulted from values of k1 and k3 corresponding with values of k1 and k3 that were found to yield chaos in figure 6. Similarly, periodic attractors resulted from values of k1 and k3 that were found to yield periodicity in figure 6.
Bifurcations of the concentration of reactants over k1 (see figure 3) or k3(see figure 4) both illustrate a period-doubling route to chaos. This property appears from either decreasing k1 or increasing k3. The period-doubling behavior on k1 appears independently of the choice of k3, as does the behavior illustrated in the bifurcation on k3 for any choice of k1. Furthermore, experimentation has verified the predicted behavior of the bifurcations (4, 5). Although the bifurcation diagrams display areas indicative of chaos, further testing is required to determine whether the disordered regions do indeed give rise to chaos. Sensitive dependence on initial conditions was found from calculating the timeseries data from the Olsen model differential quations. Timeseries data was calculated two times for initial conditions differing in [NADH] concentration by 1×10-10. Then the natural log of the magnitude of the difference of [O2] concentration was calculated, and a line was fit to the section with positive overall slope, if it existed (see figure 5). The slope of the regression line is therefore the Lyapunov exponent corresponding to the orbit of the oxygen concentration in the system. Tests for sensitive dependence were carried out by monitoring [O2], [NADH], [NAD.], and [Co III], while varying the difference in initial conditions among all four variables. 
Also, different values of k1 and k3 were used. The Lyapunov exponents resulting from the sensitive dependence tests verify the data in the bifurcation diagrams: no sensitive dependence (slope < 0) was found where the bifurcation diagram shows periodic behavior of [O2], while sensitive dependence was found (slope > 0) where the bifurcation diagrams predict chaos. The cutoff on positive Lyapunov exponents was taken to be 1e-4. Testing for sensitive dependence in the above manner proved to be very tedious, and is impractical for understanding how k1 and k3interact to either give rise to a steady state, periodic oscillations, or chaos. In order to characterize these interactions, a plot over k1 and k3 was constructed, displaying Lyapunov exponents for the different k values (figure 6). The black area in figure 6 represents values of k1 and k3 for which the concentration of oxygen is in a steady state. Note that when k3is increased just past the steady state area, the behavior of oxygen concentration becomes chaotic (light grey/white). Then as k1 decreases and k3increases, there is a large area of periodic behavior (grey), eventually getting to a region of chaos denoted by white and light shades of grey in the figure. As k1 continues to decrease and k3 continues to increase, there is once again a large area of well behaved, periodic oscillations in [O2]. Figure 6 supports the general trend seen in the bifurcation diagrams (figures 3 and 4): decreasing k1 leads from periodic to chaotic to periodic behavior, as does increasing k3. Furthermore figure 6 illustrates that this trend is relatively universal across values of k1 and k3. Lyapunov exponents were calculated as the slope of a linear regression of the natural log of the difference in [O2] (as described above for the sensitive dependence test). The Lyapunov exponents are a conservative approximation: the linear regression was taken from time 50 to time 250 in all cases, so in a case of extreme sensitivity to initial conditions, the Lyapunov exponent in figure 6 is too small (but still larger than 0.005 and therefore indicative of chaos). Finally, in an attempt to verify the chaotic nature of the imeseries data, a one-dimensional map was constructed (figure 7). 
The map constructed by plotting the preceding amplitude of the [O2] timeseries data against the following amplitude. The map appears to be a fractal since zooming in on an area of the map yields an ordered structure. Furthermore, in preliminary calculations (due to time constraints) with mediocre resolution, the box counting dimension of the fractal appeared to be around 1.1 or 1.2. This indicates that indeed the map is a fractal. The fact that the time-delay embedding of the successive peaks in the [O2] time series is a fractal supports the claim that the PO reaction indeed exhibits true chaos. If the orbit of [O2] was merely periodic with a very long time scale, the fractal dimension of figure 7 would be zero (just a finite collection of points). Showing that the fractal dimension of the map is greater than zero would provide strong evidence for [O2] exhibiting chaos.
Conclusions
There is significant and conclusive evidence that the PO reaction displays true chaotic behavior. This statement holds true only for certain choices of k1 and k3; however the existence of chaos in molecular reactions occurring inside an organism is a significant result. We have shown that the attractors are bounded for the Olsen model. Furthermore, for a given chaotic orbit we have seen that there exists a corresponding positive Lyapunov exponent. These two facts by themselves justify calling the orbit chaotic (this is the definition of chaos). Additionally we have shown that there indeed exists sensitive dependence on initial conditions, and that a one-dimensional map illustrating the successive change in maxima of the timeseries data for a reactant is a fractal. Finally, there is a period-doubling route to chaos. Thus it is clear that the PO reaction can exhibit chaos. One issue requiring further attention is exactly how accurate the Olsen model is at describing the chemical concentrations in the PO reaction. 
Experimentation has proven that the model describes periodic oscillations correctly: the number of different peaks in the experimental timeseries data corresponds with the predictions from the bifurcation diagrams. Similarly, experimental results mimic the chaotic behavior where predicted on the bifurcation diagrams. However since the Olsen model displays sensitive dependence for certain choices of k1 and k3 (presumably the PO reaction does also), it is impractical to expect the model to properly predict the chaotic timeseries data for an experimental procedure. Indeed that is the nature of chaos: it is extremely difficult to predict. The most important tests of how good the Olsen model is are those that test how well the model predicts periodic behavior, and how well the model predicts when chaos will arise.
Acknowledgements
I would like to thank Professor Alexander Barnett for all of his mathematical guidance. Also I thank Professor Robert Ditchfield and Charlie Ciambra for their help with the chemistry.
References
1. K. T. Alligood, T. D. Sauer, J. A. Yorke. Chaos: An Introduction to Dynamical Systems. (Springer-Verlag, New York, 1996).
2. I. R. Epstein, K. Showalter, J. Phys. Chem. 100, 13134 (1996).
3. I. R. Epstein, Physica D. 7, 47 (1983).
4. C. G. Steinmetz, T. Geest, R. Larter, J. Phys. Chem. 97, 5649 (1993).
5. L.F. Olsen, Physics Letters A. 94, 454 (1983).
6. I. R. Epstein, J. A. Pojman, An Introduction to Nonlinear Chemical Dynamics: Oscillations, Waves, Patterns, and Chaos. (Oxford University Press, New York, 1998).
7. L. Gyorgyi, R. J. Field, J. Phys. Chem. 95, 6594 (1991).
8. R. J. Field, L. Gyorgyi, eds. Chaos in Chemistry and Biochemistry. (World Scientific Publishing Company, 1993).