681
J. Phys. Chem. 1995,99, 681-686
Chaos Control in an Enzymatic Reaction A. Lekebusch, A. Fiirster, and F. W. Schneider" Institute of Physical Chemistry, University of Wiirzburg, Marcusstrasse 9-1I , 0-97070 Wiirzburg, Germany Received: August 8, 1994; In Final Form: October 13, 1994@
We apply the continuous delayed feedback method of Pyragas to the experimental control of chaos in the peroxidase-oxidase (PO) reaction. Unstable periodic 1' and l2orbits embedded in the chaotic PO attractor were stabilized in CSTR experiments. The stabilization is demonstrated by a minimum in the experimental dispersion function and by the equality of the delay time z and the period of the stabilized attractor. In model calculations we compare the Pyragas method with the discontinuous control method by Ott, Grebogi, and Yorke (OGY). Experimental noise reduces the control efficiency of the OGY method much more severely than that of the Pyragas method precluding the experimental stabilization of unstable periodic orbits by the OGY method in the present PO system.
Introduction The control of chemical chaos has recently received increased Ott, Grebogi, and Yorke (OGY)' formulated a discontinuous method of stabilizing unstable periodic orbits (UPO). Small single perturbations are imposed on the system when it is close to a fixed point in the Poincare section. The fixed point represents an unstable periodic orbit.3 In work by Petrov et aL2 a reduced map-based algorithm has been used which is appropriate for chemical systems showing low dimensional c h a o ~ . ~The . ~ system can be shifted toward the fixed point A, by applying a single perturbation per orbit to the bifurcation parameter p. The perturbations Ap are calculated according to the linear relationship
Figure 1. Diagram of continuous time delayed feedback control with a delay line D and a delay time t after ref 3. The output is the observed speciesy(t) which is equal to the concentration of 0 2 in the experiments and A in the calculations of the Olsen model.1° In the experiments the oxygen inflow rate is kf(02) = kfo(o2) F(t,z). In the Olsen model the flow rate of A is kf@) = kfo(A) F(t,t).
+
+
02input
where g is a constant and An is a point close to A,. The OGY algorithm has been used to control chaos in a number of physical systems6v7and in the Belousov-Zhabotinsky (BZ) reaction.2 Continuous control of chaos by time-delayed feedback has been formulated by Pyragas8 and carried out in the BZ rea~tion.~ At each data point in the time series, the time-delayed feedback function applies a small-amplitude perturbation to the total inflow rate or to the flow rate of a single species. These perturbations are calculated as the difference between the previous (time-delayed) signal y(t-z) and the actual signal y(t)
F(t,z) = K[jt(t-z) - y(t)] = KD(t,z)
(2)
where z is the delay time, K is the feedback strength and D(t,z) is the difference between the two signals (Figure 1). The dispersion (variance) D2 represents the time average value of D2(t,z)at a constant delay time. In the numerical integrations F(t,z) is used as an additive term in the differential rate equations. When the delay time z is equal to the period of an unstable periodic orbit, F(t,z) tends to zero and only small perturbations are necessary to stabilize this UPO. Recently, unstable periodic 1' and l2 orbits embedded in the chaotic attractor were stabilized by this method in the BZ rea~tion.~ The Pyragas method has also been applied to calculations of a model (Aguda-Larter model (AL)) for the peroxidase-oxidase (PO) rea~tion.~ In addition to the stabilized orbits in the strange
* Author to whom correspondence should be addressed. @
Abstract published in Advance ACS Abstracts, December 15, 1994.
0022-3654/95/2099-0681$09.00/0
-1 02Aiquid output
liquid input 2=
Figure 2. CSTR with two inlets for the reactant solutions, gas inlet, gasfiquid outlet, and 0 2 Clark electrode. Quartz glass sides are parallel to the paper plane. Plexiglas sides and UV-vis light beam are perpendicular to the plane.
attractor, new periodic, chaotic, and steady states were produced when F(t,z) attains large values for arbitrary choices of z in the AL model. In this work we experimentally apply the Pyragas method of continuous self-controlling delayed feedback to the enzymatic PO reaction. We also compare the Pyragas method with the OGY method in calculations using the Olsen modello of the PO reaction. The experimental chaotic states of the PO reaction are prepared in a CSTR (continuous flow stirred tank reactor), where 0 1995 American Chemical Society
Lekebusch et al.
682 J. Phys. Chem., Vol. 99, No. 2, I995
a
horseradish peroxidase (HRP) catalyses the aerobic oxidation of reduced nicotinamide adenine dinucleotide (NADH) according to
P1
2NADH 11
10
1.20
12 13 0 content (%)
0.98
0.so
(3)
NAD+, also known as coenzyme 1, is B-nicotinamide adenine dinucleotide. A chemical recycling of NAD+ produced in reaction 3 may be achieved by another enzyme reaction (4), which ultimately permits the use of lower flow rates in the CSTR: NAD'
0.40 I
0.60
+ G6P [DHl_ 6PGL + NADH + H+
(4)
HB
k7
=Chaos
15
14
[Wl + 0, 2NAD' + 2H20
Figure 3. Bifurcation diagrams: (a) experimentalbifurcation diagram at pH = 5.80 as a function of the 0 2 content in the gas mixture. Higher periodic patterns (14, 15, ...) are summarized in a 1" region; (b) bifurcation diagram calculated with the Olsen model as a function of the inflow rate of A, k7. HB = supercritical Hopf bifurcation point. The chaotic state used for the control is located at k7 = 0.98.
where DH is glucose-6-phosphate dehydrogenase. Chaos in the PO reaction has been found in 1977 by Olsen and Degn in a fill reactor.12*13These chaotic oscillations were reproduced under similar conditions and a period-doubling route involvine " P1. " P2. , P4. chaos, and P3 was found.14-15 Recently,
a 2.7
9 0
C
2.4
a
R
0
2.1
1.8 0
1000
2000
3000
Time [SI
4000
5000
6000
b
d
1
0.5 n
s2
W
0
R
-0.5
-1
I
I
I
I
1
I
0
1000
2000
3000 Time [SI
4000
5000
6000
Figure 4. (a) Experimental time series of the PO reaction showing chaotic oscillations and a stabilized 1' pattern after the start of the Pyragas feedback control (at t = 2000 s) for a delay time of z = 143 s and a coupling strength of K = 0.3 mL2 mol-' s-'; (b) dispersion curve of the above time series; (c) SVD-reconstructed attractor of the chaotic time series; (d) attractor of the stabilized 1' time series.
J. Phys. Chem., Vol. 99, No. 2, 1995 683
Chaos Control in an Enzymatic Reaction TABLE 1: Olsen Model with Rate Constants1o = 0.347 kz = 250 k3=0.035 h = 20 ks = 5.35 ks = 1 x 10-5 k7 = bifurcation parameter; k-7 = 0.1 ks = 0.825
(1.1) (1.2) (1.3) (1.4) (1.5) (1.6)
B+X42X 2X42Y A+B+Y43X X-P Y-Q [XIo-X (1.7) [Y]o*A (1.8) [Blo4B
kl
a period-doubling route to chaos in a Farey sequence of the PO reaction was also observed.16 The latter chaotic state is used in the present CSTR experiments.
Experimental Section Materials. Peroxidase from horseradish (HRP,RZ 3.05, activity 270 units per mg of solid) was purchased from Sigma as a salt-free powder. Glucose-6-phosphate dehydrogenase (DH) from leuconostic mesenteroides (540 units /?-NAD+ per mg of protein) was purchased from Sigma as a biuret suspension. /?-NAD+from yeast (99%)and glucose-6-phosphate disodium salt (G6P, 99%) were purchased from Sigma. 2,CDichloropheno1 (DCP) and methylene blue (MB) were purchased from Aldrich. Reactor. A Plexiglas reactor is used with two quartz glass walls and an oxygen electrode fitted on one Plexiglas side (Figure 2). The reactor has a square base of 18 mm width and 18 mm height with an effective liquid volume of 4.3 mL. The gas volume above the liquid volume is 2.0 mL. We used an unsymmetxic magnetic stirrer of about 1.5 mL self-volume driven by a small motor with a stirring rate of 750 rpm. This was found to be a good combination of efficient mixing of the
a
I
I
1000
2000
I
I
3000
4000
liquid and a constant oxygen-transfer rate through the liquid surface. The oxygen-transfer rate was determined to be 1.11 x 10-3 s-1. Preparation of M o w Components. All reactants were dissolved in aqueous 0.1 m phosphate buffer (pH 5.8) containing 1 pmol/L MF3 and 50 pmol/L DCP. AU solutions were prepared in an argon atmosphere before each measurement. We used two gas-tight syringes (Hamilton) filled with NAD+/G6P solution and DWHRP solution with the following concentrations: syringe 1, [HRP]= 150 units/&, [DH] = 5.0 units/ mL; syringe 2, [NAD+] = 3.0 mmol/L, [MP] = 50 mmol/L. All experiments were done at pH = 5.80. Continuous Flow Conditions. A self-designed high-precision syringe pump driven by a computer-linked stepping motor was employed to control the flow rate of the reactants into the CSTR. The stepping frequency was 9 Hz producing a constant flow rate of the inflow of the two solutions at kdliquid) = 1.51 x min-' (residence time z = 67 min). The flow rate of the gaseous mixture O D 2 was constant at kdgas) = 1.1 mWs in all experiments. The 0 2 content in the O D 2 mixture was used as the bifurcation parameter. Temperature. The syringes containing the input chemicals, the tubing, and the reactor were thermostated at 25 "C for all experiments. Detection. The time series of NADH absorption (360 nm detection wavelength) and compound ID (418 nm) were mesured with a diode array spectrophotometer (Hewlett Packard 8452A). The time series of the oxygen electrode potential were measured with an oxygen-selective Clark electrode and a microprocessor oximeter (WTW Oxi 96). Thus the NADH concentration, the compound III concentration, and the oxygen concentration are
I
2.3
i
2
0
8
1.7
1.4
0
Time [SI
5000
6OOO
1 -
0.5
-
h
5
v
0 -
r 2 -0.5
-
-1
0
I
,
1000
2000
I
I
I
3000
4000
5000
Time [SI
I
I
I
m
Figure 5. (a) Experimental time series of the PO reaction showing chaotic oscillations and a stabilized l 2 pattem after the start of the Pyragas feedback control (at t = 2000 s) for a delay time oft = 200 s and a coupling strength of K = 0.3 mLz mol-' s-l; (b) dispersion curve of the above time series; (c) SVD-reconstructed attractor of the chaotic time series; (d) attractor of the stabilized l2 time series; (e) part of the chaotic attractor (550-750 s) representing the l2 UPO; (f) part of the stabilized attractor (3400-3600 s).
Lekebusch et al.
684 J. Phys. Chem., Vol. 99, No. 2, 1995
a 2.9
22
1
I
t
I
2o
I
I
t
I
I
I
I
44°4A
N
a
I 0
I
I
I
I
I
1
1000
2000
3000 Time [SI
4000
5000
6ooo
0
I
I
I
I
I
I
I
20
4 I
I
6
30 40 50 Delay-Time [time units]
60
70
Figure 8. Dispersion Dz versus delay time t for the Olsen model. The coupling strength is K = 0.015 and the rate constant is k7 = 0.98. All other parameters are given in Table 1. For simplicity some calculated bistabilities are not shown.
b I
, 10
I
I
0
500
I
I
I
8
6 3
2
4 2
1000
0
2000
3000 Time [SI
4000
5000
6ooo
Figure 6. (a) Experimental time series of the PO reaction showing chaotic oscillations and a generated l 2 pattem after the start of the Pyragas feedback control (at t = 2000 s) for a delay time of t = 143 s and a coupling strength of K = 0.5 mL2mol-' s-I. The period of the newly generated orbit is T = 200 s; (b) dispersion curve of the above time series.
0
1000
1500
2000
2500
Time units
Figure 9. (a) Calculated time series of the Olsen model showing chaotic oscillations and a stabilized 1' pattem after the start of the Pyragas feedback control (at t = 500 time units) for a delay time of t = 29 time units and a coupling strength of K = 0.015;
Pyragas experiments y(t) corresponds to the 0 2 concentration in the reactor. The feedback function F(t,t) acts specifically on the inflow of gaseous 0 2 .
, 0
0.2
u
, 0.4
I
I
0.6
0.8
2 Coupling Strength K [mL mol-'s-']
1
Figure 7. Dispersion Dz versus coupling strength K for a delay time of t = 200 s (I2 pattem) in the PO experiment.
detected simultaneously with a sampling rate of 2 Hz as shown e1~ewhere.l~ The figures show the NADH absorption time series only. Oxygen Flow and Control. Two computer-linked mass flow controllers (MKS Type 1259C) for 0 2 and N2, respectively,were employed to mix and regulate any chosen gas flow rate and 02 content. The imposed feedback functions were calculated by an AT-compatible computer. Figure 3a shows the experimental bifurcation diagram of the PO reaction using the above parameters. With increasing oxygen flow, a period-doubling route to chaos in a Farey sequence occurs.16 For each control experiment a bifurcation diagram was first constructed. A relatively wide range of chaos was observed. In the present
Model Calculations The Olsen modello (Table 1) is used in the present calculations, because it shows a similar scenario (a Farey sequence including chaos through period doubling) as in our experiments (Figure 3b). The Olsen model is a dynamical four-variable model that consists of eight steps. The variables are A = 02 (corresponding to y(t)),B = NADH, and X and Yare the radical intermediates NAD' and compound 111, respectively. All rate constants and starting concentrations of the employed chaotic state are given in Table 1 (bifurcation parameter k7,, = 0.98). All integrations were done using the Gear alg0rithm.'~7'~The controlled inflow parameter is the rate constant k7, which represents the transport of gaseous oxygen into the reaction solution. The observed species is the variable A. During control with the Pyragas method F(t,t) is calculated according to eq 2. The oxygen inflow is the sum h = F(t,t), where k70 is the rate constant corresponding to the free-running chaos. For the implementation of the OGY algorithm, a one-dimensional map is constructed from the chaotic time series of the oxygen concentration A. A, and A,, are the fixed point and the actual point on the one-dimensional map, respectively. The perturbation Ap is calculated according to eq 1 and added to k70. The perturbation is applied for 5 s on each retum.
+
Results and Discussion Pyragas method. Figure 4a shows the chaotic time series (13% 0 2 in the O n 2 mixture) of the free-running PO reaction
J. Phys. Chem., Vol. 99,No. 2, 1995 685
Chaos Control in an Enzymatic Reaction
normalized dimensions (X,Y,Z axes). Next we show a l2time series (Figure 5a) which has been stabilized with z = 200 s and K = 0.3 mL2 mol-' s-l. It is part of the same chaotic state as shown in Figure 4c. The Pyragas feedback algorithm is started at 2000 s and after a 800 s transience time a l2 state with a period of 200 s is stabilized. The corresponding dispersion curve (Figure 5b) and attractors (Figure 5c,d) are also shown. The unstable l2 orbit (UPO) of the chaotic attractor is shown in Figure 5e and its stabilized l2orbit is displayed in Figure 5f. It is seen that the stabilized l2attractor is similar to the U p 0 (12). Therefore it is concluded that the orbit is indeed a stabilized periodic orbit and not a newly generated periodic state. During control the dispersion D2 decreases from 0.246 (chaos) to 0.011 (12). If an improper value for K had been chosen, newly created periodic states would have been obtained. This is shown in Figure 6a which represents a generated l2state (K = 0.5 mL2 mol-' s-'). The chosen delay time was 143 s and the coupling strength has been increased from K = 0.3 to 0.5 mL2mol-' s-l. The feedback algorithm was started at 2000 s and after a transience time of 400 s l2 oscillations were obtained. The dispersion (variance) increased from 0.082 (chaos) to 0.320 (12) during control (Figure 6b). Furthermore, this l2 attractor (not shown) is not embedded in the chaotic attractor in Figure 5c. The range of possible K values for the stabilization of a UP0 is very narrow in the experiment, as exemplified by the minimum in Figure 7. The reason for the stabilization of 1' and l2orbits (and not a P1) orbit in the experimental Pyragas method may be the higher population of these unstable orbits than that of the unstable P1 orbit in the strange attractor as seen from the experimental time series (Figure 4a). Applying eq 2 to the Olsen model we observed stabilized 1 and (1')2 UPOs at z = 29 time units and 55 time units,
N
n
0.005
0
0.01
0.015 0.02 0.025 Coupling Strength K
0.035
0.03
Figure 10. Dispersion Dzversus coupling strength K for a delay time oft = 29 time units (1' pattern) in the Olsen model.
with the oxygen concentration as the controlling species. The control algorithm is initiated at 2000 s, with a delay time z of 143 s and a coupling strength of K = 0.3 mL2 mol-' s-'. z has been obtained from one of the major peaks of the chaotic Fourier spectrum (not shown), which is close to the frequency of the free-runuing 1 oscillations. Following a 100 s transience time, a stabilized 1' pattern with a period of 140 s is obtained. After switching off the control algorithm at 4000 s the system immediately returns to its previous chaotic state. Figure 4b shows the corresponding difference curve D(t,z). During control the time-averaged dispersion D2 decreases from 0.078 (chaos) to 0.010 (1'). The experimental chaotic attractor (Figure 4c) was reconstructed by the SVD method.20 It is seen that the stabilized 1' attractor (Figure 4d) is embedded in the chaotic attractor. The SVD attractor reconstruction is used since it filters instrumental (noninteractive) noise by an orthonormalization algorithm.21 The attractors are plotted in the first three d 8
6 3
7
4
2
0
500
1000 Tune units
1500
2000
b
0.97
I
0
I
500
I
lo00 Time units
I
1500
2000
Figure 11. (a) Calculated time series of the Olsen model showing chaotic oscillations and a stabilized P1 pattern after starting the OGY control algorithm (r = 500 time units); (b) plot of flow rate perturbation pulses for the above time series; (c) attractor of the stabilzied P1 time series; (d) attractor of the chaotic time series.
686 J. Phys. Chem., Vol. 99, No. 2, 1995
Lekebusch et al.
respectively, as seen by the corresponding minima in the D2 plot (Figure 8). Newly created periodic and aperiodic states were also observed. Figure 9a shows the chaotic time series of the free-running Olsen model together with the Pyragas control which is started at 500 time units with a delay time of z = 29 time units and a coupling strength of K = 0.015. There is qualitative agreement between experiment and model. The range of possible values for K is very narrow in the calculations (Figure 10) and in the experiments. OGY Method. With the OGY method the experimental stabilization of P1, 1' or l2 orbits was not possible in the PO reaction. We offer two reasons: The P1 U P 0 is only sparsely populated in the experimental chaotic state (see above). Furthermore, the 1' and l2orbits require corresponding onedimensional maps (n versus n 2 and n 3, respectively) which turned out to be very noisy, so that a precise intersection of the map with the bisectrix could not be achieved. More importantly, the OGY method applies only a single control process per orbit when the system is close to the fixed point. Experimental noise will tend to reduce the effect of a single control event. On the other hand, the Pyragas method applies a large number of control events. In fact, a control process is applied at each single data point according to eq 2. Thus several hundred control events occur per orbit in the experiments. Therefore, the Pyragas method is considerably less sensitive towards experimental noise than the OGY method as verified by model calculations.22 In model calculations of the Olsen model a stabilization of a P1 pattern is achieved using the OGY algorithm (Figure l l a ) in the absence of any superimposed noise. The control algorithm is started at 500 time units. After a transience time of about 100 time units a P1 pattern is stabilized. The control algorithm is switched off at 1500 time units. The original chaotic state returns after about 60 time units. Figure 1l b shows a plot of the flow rate perturbation pulses versus time. The amplitude of the pulses decreases rapidly and after 100 time units small perturbations of only about 0.3% of k7 are necessary to keep the system on the stabilized orbit. The stabilized P1 attractor (Figure 1IC)fits nicely into the chaotic attractor (Figure 1Id).
+
+
Conclusion In the experiment unstable 1' and l2 orbits embedded in a strange attractor of the PO reaction have been stabilized by the self-controlling feedback algorithm of Pyragas. This was demonstrated by a minimum in the dispersion function and the equality of the delay time z with the period Tsof the stabilized oscillations. The Pyragas algorithm is very sensitive to the choice of the delay time z and the coupling strength K. Stabilizationcan be achieved only in a narrow range of K values. With all other K values new orbits are generated. The experimental results are compared with calculationsof the Olsen
model which are found to reproduce the experimental results with qualitative agreement. The more realistic Aguda-Larter model also shows generic similarities with these experiment^.^ With a discontinuous algorithm based on a method by Ott, Grebogi, and Yorke, the experimental stabilizationof UPOs was not successful in the PO reaction due to the great sensitivity of the OGY method toward the existing experimental noise. However, in the Olsen model a stabilization of P1 by the OGY method could readily be achieved in the absence of any imposed interactive noise. In agreement with earlier calculation^,^^ the speed of convergence to the stabilized orbit is limited by the delay time in the linear Pyragas method. The Pyragas control shows a longer transience time than the OGY method. This transience time can be shortened by the application of a nonlinear neural net approach to the OGY method.23
Acknowledgment. We thank the Volkswagen Stihng, the Deutsche Forschungsgemeinschaft and the Fonds der Chemischen Industrie for partial support of this work. References and Notes (1) Ott, E.; Grebogi, C.; Yorke, J. A. Phys. Rev. Lett. 1990,64,1196. (2) Petrov, V.; Gaspar, V.; Masere, J.; Showalter, K. Nature 1993,361, 240. (3) Romeiras, F. J.; Grebogi, C.; Ott,E.; Dayawansa, W. P. Physica 1992,58D, 165. (4) Peng, B.; Petrov, V.; Showalter, K. J. J . Phys. Chem. 1991,95, 4957. (5) Petrov, V.;Peng, B.; Showalter, K. J. J. Chem. Phys. 1992,96, 7506. (6) Roy, R.; Murphy, T. W.; Maier, T. D.; Gills, Z.; Hunt, E. R. Phys. Rev. Lett. 1992,68, 1259. (7) Hunt, E. R. Phys. Rev. Len. 1991,67, 1953. (8) Pyragas. K. Phys. Len. 1992,A170,421. (9) Schneider, F.W.; Blittersdorf, R.;Farster, A.; Hauck, T.; Lebender, D.; Muller, J. J. Phys. Chem. 1993, 97, 12244. (10) Olsen, L. F. Phys. Lett. 1983,A94, 454. (11) Yamazaki, I.; Yokota, K. In Biological and Biochemical Oscillations; Chance, B., Pye, E. K., Gosh, A. K., Hess, B., Eds.; Academic Press: New York, 1973; p 109. (12) Olsen, L. F.;Degn, H. Nature 1977,217, 1047. (13) Olsen, L. F.;Degn, H. Biochim. Biophys. Acta 1978,527,212. (14) Geest, T.; Steinmetz, C. G.; Larter, R.; Olsen, L. F.J . Phys. Chem. 1992,96, 5678. (15) Larter,R.; Olsen, L. F.; Steinmetz, C. G.; Geest, T. In Chaos in Chemisrry and BiochemistTy; Field, R. J., GySrgyi, L., Eds.; World Scientific: New Jersey, 1993;p 175. (16) Hauck, T.; Schneider, F. W. J. Phys. Chem. 1994,98,2072. (17) Hauck, T.; Schneider, F. W. J. Phys. Chem. 1993,97,391. (18) Gear, C. W. Numerical Initial Value Problems in Ordinary Differential Equurions; Prentice-Hall: Englewood Cliffs, NJ, 1971; p 209. (19) Hindmarsh, A. C. Gear: Ordinary Differential Equations System Solver, VCID 2001, rev 3, Dec 1974. (20) Broomhead, D. S.; King, G. P. Physica 1986,200, 217. (21) Schneider, F. W.; Miinster, A. F. J . Phys. Chem. 1991,95,2130. (22) Muller, J.; Lebender, D.; Schneider, F.W., manuscript in preparation. (23) Lebender, D.; Miiller, J.; Schneider, F.W. J . Phys. Chem., submitted for publication. JP9420842
