Principal Component Analysis of Dynamical Features in the

Inherent variance due to oscillations in the peroxidase- oxidase (PO) reaction was studied using principal com- ponent analysis (PCA). The substrates ...
0 downloads 0 Views 213KB Size
Download PDF
Anal. Chem. 2000, 72, 1381-1388

Articles

Principal Component Analysis of Dynamical Features in the Peroxidase-Oxidase Reaction Ewa S. Kirkor and Alexander Scheeline*

School of Chemical Sciences, University of Illinois at UrbanasChampaign, 600 South Mathews Avenue, Urbana, Illinois 61801 Marcus J. B. Hauser†

Physical Biochemistry Group, Institute for Biochemistry, Odense University, Forskerparken 10, 5230 Odense M, Denmark

Inherent variance due to oscillations in the peroxidaseoxidase (PO) reaction was studied using principal component analysis (PCA). The substrates were oxygen and reduced nicotinamide adenine dinucleotide (NADH). Horseradish peroxidase (HRP) catalyzed the reaction. The concentration of a cofactor, methylene blue (MB), was varied, and 2,4-dichlorophenol was kept constant. Increase in the NADH influx was used to change the reaction dynamics from periodic to chaotic. The reaction space was abstracted to the most significant, mutually independent, pairs of absorption and kinetic basis vectors (principal components). Typically, two significant principal components were extracted from the periodic time series and three from the chaotic data. The PCA models accounted for 70-97% of experimental variance. The greatest fraction of the total variance was accounted for in experiments exhibiting periodic dynamics and less than 25 nM MB. More MB induced an increased contribution of NADH to the PO oscillator variance, as did increased NADH influx. A simulated absorption time series, computed from a mass-action model of the chemistry, was analyzed by PCA as well. The comparison of simulation with experiment indicates that the chemical model renders the time series for HRP oxidation forms with fidelity, but incompletely represents NADH chemistry and other salient processes underlying the observed dynamics. The peroxidase-oxidase (PO) oscillator involves peroxidasecatalyzed oxidation of reduced nicotinamide adenine dinucleotide (NADH) with oxygen, aided by cofactors, methylene blue, and certain phenols and anilines. Periodic1 and chaotic2 oscillations of the horseradish peroxidase (HRP)-catalyzed example of the † Current address: Institut for Experimental Physics, Department of Biophysics, Otto von Guericke University, Universitatsplatz 2, 39106 Magdeburg, Germany. (1) Yamazaki, I.; Yokota, K.; Nakajima, R. Biochem. Biophys. Res. Commun. 1965, 21, 582-586. (2) Olsen, L. F.; Deng, H. Nature 1977, 267, 177-178.

10.1021/ac990957o CCC: $19.00 Published on Web 02/17/2000

© 2000 American Chemical Society

reaction have been a subject of much study.3 The complex mechanism of the PO reaction has been most successfully modeled as a chemical switch, alternating between two distinct chemical pathways.4 The best published model5 is semiquantitative. It allows for only partial, chemically plausible representation of processes giving rise to the observed wealth of dynamics under variation of any of several control parameters: the flow of substrates, changes in pH, or initial concentration of the catalyzing enzyme. Large reaction systems are encountered in many areas of chemistry, such as combustion, polymerization, or networks of biochemical reactions. Large sets of reacting species lead to complex reaction mechanisms, which frequently manifest nonlinear behaviors. These lead to a multiplicity of states whose stability depends on values of many parameters. In studying complex systems such as oscillatory reactions, determining the number, identity, and transformation pathways among the components presents substantial difficulties. Experimental approaches to system characterization, if carried out by factorial design, quickly grow prohibitively laborious even for systems having only a few variables.6,7 One means to classify variables in complex reaction systems is to apply small perturbations to an asymptotic dynamic state.8 There is always a chance that such methods may lead to qualitative, but unintentional, changes in system dynamics by shifting the system into an adjacent basin of attraction, unless an infinitesimally small perturbation is employed. Small perturbations usually result in small, transient responses, which may be difficult to precisely observe due to large relative measurement error. An example of the application of such a perturbative approach to the peroxidase system of interest here has been (3) Scheeline, A.; Olson, D. L.; Williksen, E. P.; Horras, G. A.; Klein, M. L.; Larter, R. Chem. Rev. 1997, 97, 739-756. (4) Olson, D. L.; Williksen, E. P.; Scheeline, A. J. Am. Chem. Soc. 1995, 117, 2-15. (5) Bronnikova, T. V.; Fed′kina, V. R.; Schaffer, W. M.; Olsen, L. F. J. Phys. Chem. 1995, 99, 9309-9312. (6) Meglen, R. R. J. Chemom. 1991, 5, 163-180. (7) Arkin, A.; Shen, P.-D.; Ross, J. Science 1997, 277, 1275-1279. (8) Chevalier, T.; Schreiber, I.; Ross, J. J. Phys. Chem. 1993, 97, 6776-6787.

Analytical Chemistry, Vol. 72, No. 7, April 1, 2000 1381

Table 1. Reactions Used in Simulating PO Oscillations, Largely on the Basis of the BFSO Model5 a rxn no.

reaction

rate const (M-1 s-1)c

r1 r2 r3 r4 r5 r6 r7 r8 r9 r10 r11 r12 r13 r14 r15

NADH + O2 f NAD+ + H2O2 Per3+ + H2O2 f CpI CpI + NADH f NAD• + CpII CpII + NADH f NAD• + Per3+ NAD• + O2 f NAD+ + O2•Per3+ + O2•- f CpIII O2•- + O2•- f O2 + H2O2 NAD• + CpIII f CpI + NAD+ NAD• + Per3+ f Per2+ + NAD+ Per2+ +O2aq f CpIII NAD• + NAD• f (NAD)2 NADH f NADHprod NADHstock f NADH O2gas f O2aq O2aq f O2gas

k1 ) 3.0 k2 ) 1.8 × 107 k3 ) 4.0 × 104 k4 ) 2.6 × 104 k5 ) 2.0 × 109 k6 ) 1.7 × 107 k7 ) 2.0 × 107 b k8 ) 1.3 × 108 k9 ) 1.8 × 106 k10 ) 1.0 × 105 k11 ) 5.6 × 107 k12 ) 5.56 × 10-5 s-1 k13 ) 1.693 × 10-6 s-1 k14 ) 1.69 × 10-4 s-1 k15 ) 6.0 × 10-3 s-1

[HRP]0 ) 1.5 µM

[NADH]stock ) 0.1 M [O2gas]eq ) 426 µM [O2aq]eq ) 12 µM

a Rate constants are specific to pH 5.1; proton balance is not shown. b An effective rate constant encompassing both O -• and HO • reacting 2 2 forms in equilibrium appropriate for pH 5.1. c Except reactions r12-r15.

published.9,10 A conclusion of the cited studies was that a minimum of three variables is required to describe the system dynamics. Is this correct, and can such a conclusion be drawn from observation of native variability rather than through introduction of external perturbation? The present paper attempts such characterization. Principal component analysis (PCA) is a method that allows for objective estimation of the minimum complexity model of a chemical reaction system11 without any initial assumptions regarding the phenomena under study.12,13 The experimental data are viewed as an abstraction of the system’s properties to a vector space. The analysis is done by decomposition of an appropriate data matrix onto matrices of linearly independent, orthogonal eigenvectors,14-18 selection of statistically significant basis vectors and inspection of their properties. The experimental data decomposition is accomplished with the aid of multivariate least-squares optimization.19 The PCA regression model leads to the representation of the data matrix X as a sum of a product of eigenvector matrices T and P′ (transposed matrix P) and the matrix of residuals E in the following form:

X ) TP′ + E

(1)

The rank of the set of principal components that accounts for system variance other than measurement method error or random noise is an indication of minimum system complexity. If correlations exist among the contributing species, the principal components represent linear combinations of their signals. The number of principal components may reach the number of species in the reaction system only if all the species are mutually independent and contribute detectable, nondegenerate signals to analyzed data. Ultimately, measurement noise limits the number of discernible reaction space structures and their uniqueness.18 In this study, we analyzed time series of absorption spectra of the PO oscillator. PCA-generated matrices T and P′ correspond to absorption and its temporal amplitude basis vectors. The absorbance is proportional to the concentration of species according to the Lambert-Beer law. Interpretation of these principal components requires independently derived chemical, spectral, 1382

Analytical Chemistry, Vol. 72, No. 7, April 1, 2000

and kinetic data,20 only some of which are available.3,21,22 The PO oscillator contains species that do not absorb significantly within the observed spectral window. These are the modifiers, methylene blue and 2,4-dichlorophenol (DCP), as well as oxygen-derived species, HO2•, O2•-, and H2O2. With regard to these “silent” species, our data are underdetermined. For comparison with analysis of experimental data, a PCA model of simulated, periodic, absorption time series was constructed on the basis of a chemically realistic model of the PO oscillator. The search for the most experiment-like, simulated, absorbance time series encompassed the most recent models of Olson et al.,4 Hung et al.,10 Bronnikova et al.,5 and Larter and Hemkin.23 The simulated data represent periodic dynamics calculated with the best measured elementary reaction rate constants available.3 A fitted rate constant was employed only for reaction r8 (Table 1), where no experimental data were available. DCP was not included in the reaction scheme since its role in the reaction’s mechanism has not yet been definitely established. Recent studies have shown that numerous phenolic compounds are able to induce complex dynamics in the PO system.22,24 The degree of dynamical complexity is governed by the redox potential (9) Hung, Y. F.; Ross, J. J. Phys. Chem. 1995, 99, 1974-1980. (10) Hung, Y. F.; Schreiber, I.; Ross, J. J. Phys. Chem. 1995, 99, 1980-1987. (11) Vajda, S.; Valko, P.; Turanyi, T. Int. J. Chem. Kinet. 1985, 17, 55-81. (12) Martens, H.; Naes, T. Multivariate Calibration, 2 ed.; John Wiley & Sons: Chichester, New York, Brisbane, Toronto, Singapore, 1991. (13) Joliffe, I. T. Principal Component Analysis; Springer-Verlag: New York, 1986. (14) Vajda, S. Prace Naukowe Instytutu Cybernetyki Technicznej, Politechniki Wroclawskiej.; Wroclaw Technical Univ, Inst of Engineering Cybernetics, Bierutowice: Wroclaw, Poland, 1985; pp 227-254. (15) Vajda, S. IEEE Trans. Autom. Control 1987, 2, 182-184. (16) Vajda, S.; Valko, P.; Godfrey, K. R. Automatica 1987, 23, 707-718. (17) Vajda, S.; Rabitz, H. IEEE Trans. Autom. Control 1989, 34, 220-223. (18) Vajda, S.; Rabitz, H. J. Phys. Chem. 1994, 98, 5265-5271. (19) Martens, H.; Naes, T. Multivariate Calibration; John Wiley & Sons: Chichester, New York, Brisbane, Toronto, Singapore, 1991; pp 97-101. (20) Sans, D.; Nomen, R.; Sempere, J. Comput. Chem. Eng. 1997, 21, S631S636. (21) Klemm, A.; Steiner, T.; Flotgen, U.; Cumme, G. A.; Horn, A. Methods Enzymol. 1997, 280, 171-180. (22) Hauser, M. J. B.; Olsen, L. F. Biochemistry 1998, 37, 2458-2469. (23) Larter, R.; Hemkin, S. J. Phys. Chem. 1996, 100, 18924-18930. (24) Kummer, U.; Hauser, M. J. B.; Wegmann, K.; Olsen, L. F.; Baier, G. J. Am. Chem. Soc. 1997, 119, 2084-2087.

for the formation of the respective phenoxyl radicals.22 Mass spectrometric investigations revealed that the phenolic compounds are not consumed in the course of the PO reaction. On the basis of these findings and on the reduction potentials of the enzymic and NADH species, it has recently been suggested that the phenolic compounds mediate electron transfer from NADH to CpI, CpII, and superoxide.22 Other explicit mechanistic pathways accounting for the chemistry of DCP in the PO system have been suggested in two reaction models.10,23 EXPERIMENTAL SECTION The experimental parameters for this version of the PO reaction have been published.25 Only information necessary for analysis of experimental data and construction of simulated data sets is reviewed here. The reaction occurred in a semi-batch reactor, steadily supplied with NADH and oxygen. In periodic experiments, 0.1M NADH in water was infused at 60 µL/h into a 10-mL solution of HRP (1.2-1.5 µM), DCP (25 µM), and variable concentrations of methylene blue (MB). Higher NADH influx was used to change the reaction dynamics from periodic to chaotic. The concentration of MB ranged from 0 to 200 nM. Data monitoring commenced with the initiation of oxidation. Absorbance spectra were measured from 350 to 600 nm with 1-nm resolution. Oxygen concentrations were sensed with a Clark electrode, both with frequency 0.5 Hz. Once the dynamic pattern was established, the data were collected for 3000 or 4000 s, starting no sooner than 5000 s and no later than 8000 s from initiation of oxidation. Analyzed, periodic sets start at the same phase in the oxygen cycle and contain only full cycles (at least 20 oscillations) in that time span. Labeling of Data. For compactness, the simulated data are labeled S, the periodic measurement series are labeled P, and the chaotic series are labeled C. The last two are annotated by MB concentration in nanomolar; e.g., the notation P25 designates the periodic time series obtained using 25 nM MB in the reaction mixture. The P set time series was separated into two data subsets, labeled D for “decrease” and I for “increase” for analysis with respect to oxygen concentration changes. Principal Component Analysis. PCA was conducted using The Unscrambler (version 5.03, Camo, Norway) software package. The data were represented as difference spectra from the mean spectrum for the given time series. Wavelength varied along columns and time along rows of each data matrix. Such organization leads to spectral variance in absorbance as regression variance parameter, absorbance eigenvectors appearing as constituents of the matrix T, and amplitude (kinetic) eigenvectors as constituents of the matrix P in eq 1. Cross-validation was carried out using five randomly selected subsets of each data matrix. Design saturation, identified by the first local minimum in the residual variance, indicated the number of significant principal components for specific reaction conditions. Each fitting was carried out with two weightings, one with data of equal weight and the second weighted by the inverse of the data’s standard deviation. Avoiding temporal distortion, the first case was better suited for establishing kinetics of amplitude eigenvectors. The second case yielded less noisy absorbance eigenvectors, otherwise practically identical,

aside from a wavelength-independent scaling factor. Variable weighting was applied to determination of the spectral waveforms of the absorbance eigenvectors, as it renders the analysis less random noise dependent. Since the data are represented as differences from the mean spectrum, the eigenvector waveforms are recovered as difference spectra of the species participating in the chemical pathway of a given rank. In addition, the periodic data were analyzed with respect to phase of change in oxygen concentration. Each oxygen consumption phase (D), lasting from 60 to 72 s, contained from 30 to 36 spectra. The oxygen accumulation phase (I) took from 100 to 120 s and it contained from 50 to 60 spectra. Subdivision of the P set time series yielded two data subsets with the same total variance amplitude and signal-to-noise (S/N) ratio. Limits of Regression Accuracy. The S/N ratio was established for the periodic data by comparing the envelopes and amplitudes of adjacent oscillations in a series. We assumed that the difference in amplitude and signal distribution between two adjacent, periodic, oscillations is due to measurement error (an assumption obviously inapplicable for chaotic data). Short-term deterioration of the reaction mixture is thus included with instrumental factors as noise sources. The S/N ratio for the chaotic series was extrapolated from that estimated for periodic data. the standard deviation of the measurements was presumed constant, while the total systematic absorbance variance (the signal) was proportional to the mean amplitude of oscillations, however irregular they might be. Construction of the Simulated Data Set. Upon review of mechanistic models, the BFSO model5 was selected as the basis for calculation of a simulated time series. The set of reactions and rate constants chosen is shown in Table 1. The rate constant of reaction r5 (Table 1) was taken at its experimental value26 to somewhat alleviate a problem with computed NADH accumulation. Another minor change from the BFSO model is reinstatement of the NADH degradation reaction, r12, as in the Urbanalator model.4 The rate constant of r1 corresponds to the level expected for [MB] ) 25 nM at pH ) 5.1.27 No adjustment of aciditydependent rate constants was made for the pH change from 5.1 to 6.3, the pH of the experiments considered here. The sum of concentrations of the various redox forms of peroxidase was conserved. The HRP forms are abbreviated as follows: Per3+ (ferriperoxidase), CpI (compound I), CpII (compound II), CpIII (compound III or oxyperoxidase), and Per2+ (ferroperoxidase). The integration of the rate equations for the model used the program HAVCHM,28 modified to compute in double precision under TurboPascal 6.0 (Borland International, now Inprise). The simulation diverged from observation on the extent of NADH consumption. At the end of a simulated 4000 s, the NADH concentration exceeded the average experimental value ∼5-fold. The linearly increasing excess of NADH was subtracted from the simulated concentration in the calculation of the simulated data matrix. The correction led to an average concentration of NADH approximating the experimental value. Simulated oscillations in 360-nm absorbance mimic experiment. Each species concentration-time vector was multiplied by a vector of absorbance coefficients over the experimental wavelength range. Thus, a

(25) Hauser, M. J. B.; Olsen, L. F. J. Chem. Soc., Faraday Trans. 1996, 92, 28572863.

(26) Land, E. J.; Swallow, A. J. Biochim. Biophys. Acta 1971, 234, 34-42. (27) Sevcik, P.; Dunford, B. Int. J. Chem. Kinet. 1995, 27, 995-1002. (28) Chesick, J. P. J. Chem. Educ. 1988, 65, 599-602.

Analytical Chemistry, Vol. 72, No. 7, April 1, 2000

1383

Table 2. Effect of Direction in the Oxygen Concentration Gradient in an Oscillation on Incremental Calibration Variance and Final Residual Variance for PCA Models of Corresponding Subsets in P Experiment Series varianceb P set [MB]a 0 5 10 25 50 100 200

O2 decrease R1 93.8, 95.3, 94.9, 96.0, 85.4, 90.9, 88.5,

R2 5.3, 2.5, 3.4, 3.2, 14.6, 3.8, 5.7

R3 0.3 1.2 0.7 0.3 0.1 2.2

diffc

O2 increase VE 0.6 1.0 1.0 0.5 0.9 3.1 5.8

R1 78.7, 75.5, 69.1, 66.1, 76.5, 60.1, 61.8,

R2 10.6, 10.8, 25.9, 30.4, 9.5, 26.2, 12.2

R3 1.6 3.2 0.7 0.4 1.9 1.2

VE 9.0 10.5 4.3 3.1 12.1 12.5 26.0

∆Vc1 15.1 19.8 25.8 29.1 8.9 30.8 26.7

a Nanomolar concentration of MB. bR , percent of total variance and i rank designation. VE, residual variance after extraction of all the listed principal components. c ∆Vc1, difference in the percentage of the total variance modeled by the first PC for uptake vs accumulation segments of oxygen oscillations.

matrix of absorbance vs wavelength vs time was constructed for each absorbing species (Per3+, CpI, CpIII, Per2+, NADH,29 CpII,30-32 and NAD•26,33), and all relevant matrices were added to make the simulated X matrix. The absorption coefficient of NAD• was estimated at ∼2.5 × 103 at its 405-nm maximum.34 Neither DCP nor MB absorbs significantly within the observation range. The participation of DCP has been neglected. The dynamical contribution of MB was accounted for only through modification of the rate constant of NADH oxidation by oxygen (reaction r1) according to inferences by Olson and Scheeline.35 The simulated data set was subjected to the same type of PCA as the experimental data. Because simulated data are noise-free, the matrix of residuals is insignificant. RESULTS A tabulation of cumulative variance for various rank models for each data and simulation series is shown in Table 2. Saturation ranks are listed as a function of [MB]. Discussion is based on the trends exhibited in this table. Selection of Principal Components. The estimated signalto-noise ratio exceeded 90 for the P series and 50 for the C series of experiments. Occasionally a phase reversal of one of the pair of eigenvectors is encountered, which is compensated by simultaneous phase reversal of the other eigenvector. Assessment of standard error of measurement of the total variance indicates that the magnitude of variance modeled by the second principal component is decidedly above the generally accepted quantitation limits (S/N ratio of 10 for quantitation and S/N ratio of 3 for detection), and the variance modeled by the third PC is above but near the detection limit for periodic oscillation experiments (at least when the entire time series rather than subsets of the (29) Hauser, M. J. B.; Olsen, L. F. Plant Peroxidases: Biochemistry and Physiology. IV International Symposium, Vienna, Universite de Geneve, 1996; pp 8288. (30) Chance, B. Arch. Biochem. Biophys. 1952, 41, 404-415. (31) George, P. Biochem. J. 1953, 54, 267-276. (32) Critchlow, J.; Dunford, H. B. J. Biol. Chem. 1972, 247, 3714-3723. (33) Burnett, R. W.; Underwood, A. L. Biochemistry 1968, 7, 3328-3333. (34) Swallow, A. J. Photochem. Photobiol. 1968, 7, 683-694. (35) Olson, D. L.; Scheeline, A. J. Phys. Chem. 1995, 99, 1204-1211.

1384 Analytical Chemistry, Vol. 72, No. 7, April 1, 2000

series are analyzed). For the chaotic series, the first PC is well quantitated, the second is quantifiable, and the third is at least marginally quantifiable. Of the three criteria for assessing information content of the PCA result (measurement signal-to-noise ratio, the first local minimum of the residual variance, and random structure of the residual variance), measurement signal-to-noise ratio is the most stringent. The reported variance, modeled by the higher rank principal components, does not pass the experimental error criterion although it fulfills the condition of being above the first local minimum of the residual variance. The residual variance displays marginal nonrandom temporal structure (result not shown), indicating contribution of deterministic phenomena in addition to noise. The average residual variance for the P series is 6.6%, while the C series is either 26.2% of the total after the third principal component or 20.6% at the point of exhaustive modeling. The residual variance, greater than 6% of the total for P or 12% for C sets, is in excess of three standard deviations of the average measurement error. Effect of NADH Influx on Systematic, Modeled Variance. While the oscillatory dynamics changes from periodic to chaotic as influx of NADH is varied, the distribution of the systematic variance among the principal components changes as well. In the case of the chaotic data, the average contribution of the first PC drops to 53% from 75% in the periodic case. The lowest PC1 contribution to systematic variance in the chaotic set is 37.4% of the total. At equivalent MB concentration, the contribution of PC2 to total variance in C sets (18%) relatively increases by ∼20% over the periodic data (average of 14.5% of the total). Average contribution of PC3 increases by 50% from 2.1% for P to 3.1% for C sets. Effect of MB on Systematic, Modeled Variance. For time series consisting of the complete oscillatory cycle, periodic series exhibit a nearly constant contribution from the first principal component (76-81%) until the concentration of MB reaches 50 nM. Upon further increase in MB concentration, the first principal component (PC1) contribution falls to 63% and contribution of the second principal component (PC2) increases from a level of ∼15% to ∼21% of total modeled variance. This is a 25% increase in the relative importance of PC2. The role of MB was further assessed through differences between PCA results for consumption and accumulation phases of the oxygen cycle. The residual variance for the consumption phase was much smaller (average 1.8%) than for the oxygen accumulation phase (average 11.1%). The disparity in the degree of modeling (Table 2) increased with growth in MB concentration, indicating a growing difference in the effect of MB on these two phases of the reaction. Systematic Variance in the Simulated Data. PCA decomposition of the simulated data S leads to a model saturated at two principal components, the first corresponding to 99% of the spectral variance and the second corresponding to the remaining 1%. Since the simulated data are noiseless, the residual variance is insignificant. Absorbance Eigenvectors of Experimental Data. We attribute bands present in the absorbance eigenvectors to various species known to be present in the oscillator on the basis of spectral properties of the pure compounds and the prior, direct, spectral decomposition of the oscillatory data.29 A spectral distribution of the PC1 absorbance eigenvectors, typical of all

Figure 1. Typical absorbance eigenvectors rendered by PCA modeling of the HRP oscillator. Shown are first (A) and second (B) principal components of the P, C, and S series. Experiments were run at 10 nM methylene blue.

experiments, is shown in Figure 1A (See also Figure 3 A and B). The distribution of PC2 absorption eigenvectors is shown in Figures 1B and 2, and PC3 in Figures 2 and 3. The PC1 spectral envelope is remarkably constant for all experiments and the simulation. The PC1 absorbance eigenvector exhibits a minor shoulder or maximum near 355 nm, a major broad minimum near 395 nm, a major maximum at ∼420 nm, and two more minor maxima at ∼540 and 580 nm. The 355-nm absorbance is customarily attributed to NADH. The known HRP forms, Per3+ and CpI (402 and 395 nm),36 and CpIII (peaks at 420, 540, and 580 nm), can account for the other peaks. CpIII appears at opposite phase from the band tentatively assigned to ferric or CpI forms. It is in-phase with the species absorbing at 355 nm. No combination of HRP forms at appropriate phases could account for the 355nm band. The second eigenvector (Figures 1B, 2B, and 3C,D) exhibits a major maximum at ∼360-362 nm, whose amplitude and position depend on MB concentration. It is red-shifted by ∼7 nm from the 355-nm position of the first short-wavelength shoulder or band exhibited by the first principal component. If NADH and Per3+ absorbance were in phase, then addition of the spectra would produce a red-shift of the NADH maximum whose magnitude could emulate observation depending on the proportion of HRP in the combination. The contribution of presumed NADH absorbance to the PC2 absorbance eigenvector is out of phase with Per3+, CpI, or Per2+ and in phase with CpIII. The amplitude of CpIII absorption is too small to effect a shift of the “presumed NADH” peak from 340 to 360 nm or longer wavelength. The result is that (36) Metodieva, D.; Dunford, H. B. Arch. Biochem. Biophys. 1989, 272, 245253.

Figure 2. Impact of MB concentration on higher rank absorbance eigenvectors of the periodic, P, and chaotic, C, series of experiments (The traces are arbitrarily shifted to expose shape changes, ∆ bars represent magnitude of the eigenvectors).). (A, B) Rank 2 eigenvectors of the P and C sets. (C) Rank 3 eigenvectors of the chaotic sets.

the peak in 360-370-nm range does not correlate with any of the Soret bands of HRP forms at appropriate phase and amplitude, nor does it correlate with any other detectable absorbance peak. The lack of correlation points to the involvement of a previously unaccounted species absorbing in the 355-375-nm range, likely derived but different from NADH. The combination of the 355nm shoulder and the 360-370-nm band shows significant alteration in proportions and is phase correlated to the change in rate of NADH influx (Figure 2, P vs C runs). This peak is followed by a relatively shallow minimum at 400 nm that likely corresponds to Per3+. The second maximum in the PC2 absorbance eigenvector (418-420 nm, Figures 1 and 2) corresponds to CpIII. The amplitude of this peak is sensitive to MB concentration, maximizing for [MB] ) 10 nM. As concentration of the MB is increased, CpIII intensity declines, melding into the background at 50 nM. Simultaneously the 360-nm feature intensifies and shifts to shorter wavelength. The band of CpIII is at opposite phase to a band at ∼440 nm (Per2+). In PC2, Per2+ is present in quantity only in periodic runs with 5 or 10 nM MB. When significant, the PC3 absorption eigenvector is defined mostly by contributions from bands at 420 and 440 nm (Figure Analytical Chemistry, Vol. 72, No. 7, April 1, 2000

1385

Figure 3. Differences in absorbance eigenvectors corresponding to depletion (A, C, E) and accumulation (B, D, F) phases of periodic oxygen oscillations. Rank 1 (A, B), 2 (C, D), and 3 (E, F) at MB concentration indicated (in nM). (The traces are arbitrarily shifted to expose shape changes. ∆ bars measure the magnitude of the eigenvectors).

2C). It represents only a minor contribution to the variance of the complete periodic data. It becomes more important for the oxygen consumption subset of the periodic data with low MB content (Figure 3E). For the chaotic data, it exhibits more significant weight in both absolute and relative terms then the periodic sets. Analysis of Oxygen Depletion and Accumulation Phases of the Oscillator. The oxygen consumption phase of NADH oxidation is generally better modeled than is the oxygen accumulation phase. The fraction of the total variance modeled by the first principal component is always larger during the oxygen consumption phase than during the accumulation (∆Vc1 of Table 2). Conversely, the increment of variance modeled by the second PC is conspicuously larger in the accumulation phase than in the oxygen consumption phase. The fraction of total modeled variance declines as MB concentration increases above 25 nM. The PC1 absorbance eigenvectors of the two oxygen phases are alike, but inverted (Figure 3A,B). During oxygen depletion (Figure 3C), the second rank waveform has a global maximum 1386

Analytical Chemistry, Vol. 72, No. 7, April 1, 2000

at 360 nm, a minimum at ∼405 nm, a minor maximum near 420 nm, and another minimum at ∼440 nm. The “405-nm feature” was not encountered in analysis of undivided data sets. Its origin is unknown. The PC2 of the oxygen accumulation phase (Figure 3D) presents a different picture from the oxygen depletion phase (Figure 3C). The extremum at 440 nm is absent. For all the MB concentrations, the global maxima are located at 360 nm, followed by minima at ∼400 nm and secondary maxima at 420 nm. In order of increasing MB concentration in this subset of the data, features attributable to HRP intermediates diminish, and the 360-nm maximum transforms into a shorter wavelength shoulder with higher amplitude. At 200 nM MB, the short-wavelength shoulder is the only remaining feature. The third PC eigenvector is detectable only during the oxygen consumption phase and at low MB concentration. Its main spectral features are extremes at 420 and 440 nm. Eigenvectors of the Simulated Data. The rank one absorbance eigenvector of the simulated data set (Figure 1A) is

relatively earlier in the oxygen consumption phase the more MB is present in the reaction mixture. Other phase plots represent simple limit cycles whose diameters decrease as a function of MB concentration. These results indicate that increase of MB concentration accelerates the enzyme turnover, thus further increasing oxygen consumption.

Figure 4. Comparison of periodic experiments with simulations. Experimental (A) and simulated (B) kinetic eigenvectors and oxygen traces.

remarkably similar to the first PC of the experimental data in any circumstances. The waveform of the PC2 absorbance eigenvector of the simulated set (Figure 1B) qualitatively differs from the experimental. There are two major differences between simulation and experiment at this level. The simulation places almost all of NADH chemistry in PC1. Also, the simulated data lead to misrepresentation of NADH chemistry in PC2, which is shown with incorrect phase and order of magnitude smaller intensity than in the second principal component of the experimental data. In addition, the participation of Per2+ is also inadequately represented. The kinetic eigenvectors of the simulated data differ from the experimental at all ranks (Figure 4). The predicted period is ∼50% longer than observed. The rank one vector shows steep, switchtype rises and falls in the simulation, while the experiment shows gradual, sigmoidal buildup and a “switch-off” decline. Experimental PC1 and PC2 amplitude-phase relations with oxygen kinetics are opposite to the phase relations in simulation (Figure 4). Phase Plots and Oxygen Kinetics. Phase plots of PC1 amplitude vectors against oxygen concentration (Figure 5), appearing as a tilted figure eight, exhibit two loops, one in the “+ +” quadrant and the second in the “- -” quadrant of the plot. This is a feature that has not previously been extracted from experimental data, although it was reported from numerical integration of the high (seventh)-order model of the HRP oscillator by Fed′kina et al.37 The loops undergo contraction along the oxygen coordinate while retaining amplitude along the PC1 coordinate as the concentration of MB in the mixture increases. The contraction of the loop located in the “+ +” quadrant is much more pronounced than the contraction of the second in the “- -” quadrant of the plot. The floor of the oxygen oscillations increases by a few percent, while the ceiling suddenly declines by more than half, as the MB concentration exceeds 50 nM. The minimum of the PC1 kinetic eigenvector is achieved faster and (37) Fed′kina, V. R.; Bronnikova, T. V.; Ataullakhanov, F. I. Biophysics 1992, 37, 781-789.

DISCUSSION AND CONCLUSIONS We adopted an independent assessment of the measurement error as the criterion for selection of significant basis vectors (principal components), leading to stringent but specific confidence limits. The residual variance, greater than 6% of the total, is in excess of three standard deviations of the average measurement error. Some of the PCA regression models of the absorbance data become saturated at a small number of components, leaving from 5% up to almost 30% of the systematic variance unmodeled. A nonrandom temporal structure of the residual variance calculated for these data sets (data not shown) indicates that modeling limits encountered are due to limits to the representativeness of signals with respect to the reacting species. Since the available data matrices lack signals from several species known to participate (O2, H2O2, and the modifiers among them), we restricted ourselves to qualitative analysis of the dominant features of the principal components. It is noteworthy that the number of principal components found in reducing all of our data is the same as the dimension of the embedding space for the same experiments. Periodic data are described by two principal components and can be described dynamically in a two-dimensional phase space. A chaotic dynamics is induced by increased supply of NADH to the PO reaction. With this transition, the amplitude of the Per2+ signal in the data increases. In effect, the chaotic data require the third principal component for a reasonable PCA model. Generally, chaos requires at least a three-dimensional embedding space, but may in some instances (though not with the PO reaction) require even higher dimensionality. There have been at least three prior instances in which it has been shown that chaos in the PO/NADH system occurs in a three-dimensional space.5,9,10,37 The earliest of these (Fed′kina et al.37) was based on simulation of chemically realistic models. Ross’s work9,10 was a combination of theory and experiment based on perturbation of an oscillating system. The most recent prior work, by Bronnikova et al.,5 analyzed both simulated and experimental data. All used a dynamical systems paradigm in approaching data analysis. Here, using a statistical construct dating to the 1930s, we have computed the same rank for the system dynamics without having to perturb it. The rank found from PCA can reach at most the dimension of an embedding space.38 In the absence of accidental degeneracy, the number of principal components represents the minimal rank of a vector space spanning the data, and the embedding dimension enumerates irreducible basis vectors spanning the embedding space. In such case, it is theoretically expected that the two ranks are equal.38 Methylene blue effects several aspects of the reaction system. Directly apparent is acceleration of oxygen consumption, decrease in oscillation amplitude with shorter oscillation period, and decrease in cyclical accumulation and consumption of Per2+ at (38) Gilmore, R. Rev. Mod. Phys. 1998, 70, 1455-1529.

Analytical Chemistry, Vol. 72, No. 7, April 1, 2000

1387

Figure 5. First rank kinetic eigenvectors’ phase portraits vs oxygen concentration. Median-centered principal components of the P series. (A) [MB] (nM): 0, black; 5, green; 10, red. (B) [MB] (nM): 25, blue; 50, light green; 100, magenta.

high concentrations. PCA modeling of all data matrices confirms these insights and places the effect of MB at the level of the second and third ranks both in amplitude and in relative composition. The above findings confirm the original supposition that MB regulates the contribution of Per2+ to the reaction system.29 It reveals that the decline in Per2+ participation is not monotonic; the Per2+ contribution attains a maximum at a low concentration of MB, declining as MB concentration exceeds some threshold. Clearly, MB increases consumption of NADH. If increased consumption of NADH leads to increased production of NAD•, then the decline in appearance of Per2+ might result from unavailability of Per3+. An acceleration of reduction of CpIII to CpI may become a viable short-cut in the peroxidative cycle, also explaining the decline in the feature attributable to CpIII. In parallel with acceleration of the HRP turnover via the abbreviated peroxidative cycle, one might expect increasing accumulation of (NAD)2 if, indeed, the dimer is unreactive. Certainly, it is necessary to revisit MB chemistry with NADH derivatives at the conditions of the PO reaction. The features described above have been derived by comparison of the PC absorbance vectors extracted from experimental time series with the stationary absorption spectra of chemical species known to be involved in the PO reaction. In the following, we would like to contrast these observations with data obtained from numerical simulation of a mechanistic model. This approach allows identification of the strengths and weaknesses of the commonly used mechanistic reaction models. The comparison leads us to believe that the BFSO model represents the qualitative chemistry of HRP nearly adequately. Once a correction for accumulation of NADH in the simulated reaction mixture is made, absorbance eigenvectors of the first principal component for the simulated and experimental data are indistinguishable. The (39) Bernofsky, C.; Wanda, S.-Y. C. J. Biol. Chem. 1982, 257, 6809-6817. (40) Miksic, J. R.; Brown, P. R. Biochemistry 1978, 17, 2234-2238. (41) Oppenheimer, N. In Pyridine Nucleotide Coenzymes. Part A. Coenzymes and Cofactors; Dolphin, D., Poulson, R., Avramovic, O., Eds.; John Wiley & Sons: New York, 1987; pp 324-362. (42) Klemm, A.; Steiner, T.; Cumme, G. A.; Horn, A. Anal. Biochem. 1993, 212, 375-380. (43) Avigliano, L.; Carelli, V.; Casini, A.; Finazzi-Agro, A.; Liberatore, F. Biochem. J. 1985, 226, 391-395. (44) Kirkor, E. S.; Scheeline, A., manuscript submitted.

1388

Analytical Chemistry, Vol. 72, No. 7, April 1, 2000

discrepancies appear in the temporal behavior of the amplitudes and absorbance eigenvectors of higher rank. The simulated but unobserved excess of NADH in the reaction mixture and the PCA analysis of the experiments suggest a direction for correcting the model. The spectral distribution of PC2 can be plausibly interpreted as predominantly reflecting the chemistry of NADH. These processes are responsible for ∼1 order of magnitude more systematic variance in the data than is currently allowed by the chemical model. At the level of individual species, the discrepancy may be even more serious, as NADH and cogeners have 1 order of magnitude smaller absorptivity than the HRP forms in the Soret region. The rank two absorption eigenvectors indicate significant activity of NADH-derived species. To identify this compound or compounds, several alternatives require consideration. One is β-NADH and R-NADH anomerization. Given the duration of experiments, the anomerization leads to a nonnegligible equilibration of the NADH isomers at pH 6.3. The ∼360-nm feature in PC2 is not incompatible with R-NADH absorption.21 The R-NADH anomer can react with DCP or MB as well as β-NADH can. Whether R-NADH can be utilized by HRP is an open question (though conventional wisdom is that the answer is no). HRPcatalyzed oxidation of β-NADH yields β-NAD+ without anomerization.39 Second is acid-catalyzed hydration and cyclization of NADH.40-42 Third is a NADH radical dimer, also a reducing species.33,43 We anticipate a need for investigation of NADH radical dimer reactivity in the PO oscillator. Such work is concurrently being reported.44 ACKNOWLEDGMENT Helpful discussions throughout this work with Lars Folke Olsen are deeply appreciated. This work was supported in part by NATO Grant CRG-950750. Additional funds were obtained from the U.S. National Science Foundation (Grants CHE 93-07549 and CHE 96-15739) and the Danish Natural Sciences Research Council. We thank the Carlsberg Foundation for making the diode array spectrophotometer available. Received for review August 20, 1999. Accepted January 6, 2000. AC990957O