10090
J. Phys. Chem. A 2010, 114, 10090–10096
Analysis of the Bromate-Sulfite-Ferrocyanide pH Oscillator Using the Particle Filter: Toward the Automated Modeling of Complex Chemical Systems Naoya Sato,*,† Hiroshi H. Hasegawa,†,‡ Rika Kimura,¶ Yoshihito Mori,¶ and Noriaki Okazaki§ Department of Mathematical Sciences, Ibaraki UniVersity, Bunkyo, Mito 310-8512, Japan, Center for Complex Quantum Systems, UniVersity of Texas, Austin, Texas 78712, Department of Chemistry, Ochanomizu UniVersity, Ohtsuka, Bunkyo-ku, Tokyo 112-8610, Japan, and Research Institute of Electrical Communication, Tohoku UniVersity, Katahira, Aoba-ku, Sendai 980-8577, Japan ReceiVed: May 26, 2010; ReVised Manuscript ReceiVed: August 11, 2010
This study was aimed at identifying a quantitatively accurate reaction model of the bromate-sulfilte-ferrocyanide (BSF) pH oscillator by using the simulation-based model estimation algorithm known as the particle filter. The Ra´bai-Kaminaga-Hanazaki (RKH) model proposed for the BSF system was extended by adding the protonation equilibrium of SO24 , for which the particle filter analysis was carried out to optimize the rate constants involved with reference to the measured pH oscillation data. The extended RKH model with the optimized rate constants almost completely reproduced the measured pH oscillations and the state diagram, showing the validity of the present analysis. Chemical oscillators such as the BSF system show drastic switching of the dominant reaction path, which strongly disturbs the convergence of the rate constants if the objective function is defined in a conventional manner to reflect only a single time step datum. In this study, the objective function was defined as the residual sum of squares with respect to pH taken over an interval longer than one oscillatory period, so that all of the relevant reaction steps can contribute to the objective function. This is the first report which exemplifies the effectiveness of the particle filter in the analysis of real complex chemical systems. Introduction Epstein and Pojman, in their book “An Introduction to Nonlinear Chemical Dynamics”,1 wrote “The art, and it is an art, of constructing a mechanism for a reaction complicated enough to exhibit chemical oscillations requires skill, care, and intuition.” In this paper, we try to construct a mechanism of a complex chemical reaction with the aid of information technologies based on statistical mathematics. We believe that our trial opens the door to the automated modeling of complex chemical reactions. Even for the experienced chemist, many trials and errors may occur before obtaining an acceptable model of a complex chemical system. The difficulties arise mainly for two reasons. First, it is virtually impossible to measure all of the intermediate species experimentally, which may lead to incomplete or inaccurate modeling of the reaction network structure. Second, the rate equations are not always available for all of the component reactions, which makes it difficult to give a quantitative account of the dynamical system characteristics such as period and amplitude of the chemical oscillations. In principle, unknown reaction rates should be determined experimentally from kinetic studies. However, it is not always possible to isolate a single reaction step due to the instability of intermediate reactant species. The uncertain factors such as unknown rate constants should be estimated by careful and repeated comparison of the model simulation with measured data. However, this procedure has been mostly conducted based * To whom correspondence should be addressed. E-mail: sato.naoya@ gmail.com. † Ibaraki University. ‡ University of Texas. ¶ Ochanomizu University. § Tohoku University.
on ad hoc and subjective criteria, and no systematic approach has been taken to develop the automated modeling of complex chemical systems. On the other hand, systematic approaches have been developed for the dynamic model estimation in the field of time series analysis. The classical approaches based on the state space model (SSM) are shown to be highly effective in estimating the linear system models with a Gaussian noise. For example, the Kalman filter has been used in a wide range of engineering and econometric applications.2-4 However, the SSM-based approaches are not directly applicable to the modeling of complex chemical systems due to the following reasons. First, the system equations governing the reaction kinetics are generally nonlinear. Second, unidentifiable minor reactions may contribute as a source of non-Gaussian system noise. Third, the system equations obtained from the conventional methods are in the form of recurrence equations which cannot be directly translated into the system rate law equations. A possible solution to overcome these difficulties is to employ a recently developed filtering technique called the particle filter based on the generalized state space model (GSSM). The particle filter, also known as the sequential Monte Carlo method, is the simulation-based model estimation technique proposed by Kitagawa5 and Gordon et al.6 It is widely used to estimate Bayesian models and has important applications in many fields such as econometrics, robotics, and biology.7,8 Recently, Adachi, Washio, and Motoda proposed the interesting idea of applying the particle filter to identify a coupled differential equation system representing the first-principle law of an object process.9 The performance of their method was demonstrated through the analysis of artificial data. However, their approach has not yet been applied to real systems.
10.1021/jp106700f 2010 American Chemical Society Published on Web 08/30/2010
Analysis of the BSF pH Oscillator
J. Phys. Chem. A, Vol. 114, No. 37, 2010 10091
TABLE 1: Reaction Mechanism and Rate Constants in the RKH Modela no.
reaction +
+H h
R1 R2
+ HSO3 + H h H2SO3
R3 R4 R5
2+ 3HSO3 + BrO3 f 3SO4 + Br + 3H + 23H2SO3 + BrO3 f 3SO4 + Br + 6H + 43BrO+ 6Fe(CN) + 6H f Br + 6Fe(CN) 3 6 6 + 3H2O
a
rate constant (35 °C)
rate law + k1[SO23 ][H ] (forward) k1′[HSO] (reverse) 3 + k2[HSO3 ][H ] (forward)
HSO3
SO23
k2′[H2SO3] (reverse) k3[HSO3 ][BrO3 ] k4[H2SO3][BrO3] k5[H+]/(k5′ +[H+])
k1 k1′ k2 k2′ k3 k4 k5 k5′
) ) ) ) ) ) ) )
5.0 × 1010 M-1 s-1 5.0 × 103 s-1 6.0 × 1010 M-1 s-1 1.0 × 109 s-1 5.97 × 10-2 M-1 s-1 18.0 M-1 s-1 1.5 × 10-5 M s-1 2.5 × 10-4 M
Reference 12.
In this paper, we report on the model identification of a complex chemical system through the analysis of experimental time series data by using the particle filter. As a target system, we selected the bromate-sulfite-ferrocyanide (BSF) system. The BSF system is one of the well-characterized pH oscillators discovered in the project “Systematic Design of Chemical Oscillators” lead by Epstein.1,10 The reaction mechanism of the BSF system was first proposed by Edblom et al.11 Later, Ra´bai, Kaminaga, and Hanazaki reexamined the mechanism and proposed an alternative model (RKH model), which was shown to better describe many aspects of the experimental observations.12-14 In this study, we employed the RKH model as a basic model to be optimized and extended it to precisely reproduce the time series oscillation data. To our knowledge, this is the first application of the particle filter to the model identification of real chemical systems. Methods Experimental Procedure. Reagent-grade chemicals NaBrO3, Na2SO3 (Kanto Chemical), K4Fe(CN)6, and H2SO4 (Wako Pure Chemical) were used without further purification. Distilled and deionized water was used throughout. The experiment was performed in a continuous-flow stirred tank reactor (CSTR) with a volume of 16 mL. The aqueous solutions of NaBrO3, Na2SO3, K4Fe(CN)6, and H2SO4 were prepared immediately before use and were separately fed to the reactor by using a peristaltic pump (EYELA, MP-3). The initial concentrations fixed throughout a series of experiments were [BrO3-]0 ) 75 mM, [Fe(CN)64-]0 ) 15 mM, and [H+]0 ) 15 mM, where [X]0 denotes the concentration of X in the reactor to be attained if no chemical reaction took place. The initial concentration of SO23 was varied in the range of 50 e [SO32-]0 e 160 mM. The reaction mixture was stirred vigorously by a magnetic stirrer. The normalized flow rate (inverse residence time) k0 was varied in the range of 0.025 e k0 e 0.1 min-1. The reaction temperature was kept constant at 40.0 ( 0.5 °C by circulating thermostatted water between a water jacket around the reactor and a water bath (EYELA, NTT-1200). The system pH was monitored by a calibrated glass electrode (Horiba, 6861-10C) and was recorded by a personal computer through a pH meter (HORIBA, F-23) and a 12-bit A/D converter (CONTEC, AD12-16) with a sampling rate of 2 Hz. Reaction Model. As a starting point, we employed a reaction model proposed by Ra´bai, Kaminaga, and Hanazaki (RKH)12 as given in Table 1. In this model, a pH oscillation in the BSF system is briefly explained as follows. In the presence of SO32-, reactions R1 and R3 constitute a dominant reaction channel (S channel). The sum (R1) × 3 + (R3) gives the stoichiometry 2+ 3H+ + 3SO23 + BrO3 f 3SO4 + Br + 3H
(S)
TABLE 2: Additional Reaction Mechanism and Rate Constants for the Extended RKH Model no. R6
reaction SO24
+
+H h
HSO4
rate law
rate constant (20 °C)a
+ k6[SO24 ][H ]
k6 ) 1.0 × 1011 M-1 s-1
(forward) k6′[HSO4] (reverse)
k6′ ) 1.0 × 109 s-1
a Eigen, E.; Kurtze, G.; Tamm, K. Z. Elektrochem. 1953, 57, 103–118.
which is catalytic but not autocatalytic with respect to H+. Reaction S proceeds as far as SO32- is existent, during which [H+] is kept extremely low by the rapid equilibrium (reaction R1). When SO32- is consumed, a sizable amount of H+ is liberated to generate H2SO3 through reaction R2. This triggers an autocatalytic production of H+ and hence a sharp drop of pH as given by (R2) × 3 + (R4) (D channel) 2+ 3H+ + 3HSO3 + BrO3 f 3SO4 + Br + 6H
(D) On the other hand, reaction R5 works as a negative feedback which consumes H+ steadily to bring the system back to the high-pH stage. In Table 1, the reaction rates are assumed to follow the law of mass action except for reaction R5. Reactions R3 and R4 are the composite reactions in which the first steps are considered to be rate-determining. Reaction R5 is assumed to follow the empirical rate law with a hyperbolic dependence on [H+]. In this study, we selected the rate constants k3, k4, k5, and k5′ as variable parameters to be optimized. We fixed k1, k1′, k2, and k2′ in the rapid protonation equilibria reactions R1 and R2, whose equilibrium constants are well-defined. We confirmed that small variation in the rate constants in reactions R1 and R2 does not significantly change the calculated pH waveforms. Table 2 shows an additional proton buffer system required to attain quantitative agreement with the experimental results (see the Model Identification by Particle Filter section in the Results). No measurements are available in the literature for the temperature dependence of the rate constants in reactions R1, R2, and R6. We did not attempt temperature corrections of these parameters to avoid unnecessary ambiguities. Particle Filter Algorithm. In the present analysis, a particle is defined as a set of the rate constants k3, k4, k5, and k5′. We optimized the rate constants by using the algorithm as described below. (1) A series of random particles was generated by assigning random values to the rate constants, such that each rate constant was distributed uniformly on the log-transformed scale within the range from 1/10- to 10-fold the corresponding value in the RKH model. The coupled ordinary differential equations (ODE) derived from the rate equations were numerically solved for
10092
J. Phys. Chem. A, Vol. 114, No. 37, 2010
Sato et al.
each particle. After 5000 s of reaction time, if the solution converged to an oscillatory state and the range of pH oscillations covered the pH value measured at t ) 0, the particle was selected as a candidate. We prepared M ()10 000) candidates by repeating the procedure. (2) For each particle, numerical simulation was performed for the time interval from t ) 0 to T ()2300 s). The initial concentration of H+ was set the same as that observed at t ) 0, and the other initial concentrations were appropriately determined from the numerical solution obtained in Step 1. Note that T should be set sufficiently longer than the oscillatory period (see Discussion). (3) Particle selection/reproduction was carried out by evaluating the objective (cost) function Ei, defined as the residual sum of squares (RSS) with respect to pH
Figure 1. Measured pH oscillations in the BSF system. Experimental + 4conditions: [BrO3 ]0 ) 75 mM, [Fe(CN)6 ]0 ) 15 mM, [H ]0 ) 15 -1 mM, [SO2] ) 90 mM, k ) 0.05 min , and temperature ) 40.0 °C. 3 0 0
T/τ
Ei )
∑ (pHicalcd(nτ) - pHmeasd(nτ))2
(1)
n)1
where suffix i denotes the particle number i ) 1, ..., M and τ is the data sampling interval. We used a sufficiently short sampling interval of τ ) 2 s in order for the sharp drop of pH in the autocatalytic phase to be traced with a sufficient resolution. We employed the Boltzmann distribution as a selection probability function
Pi ) C exp(-βEi)
(2)
where C is the normalization factor and β was set at 50. The M new-generation particles were reproduced according to Pi with a small Gaussian noise imposed on each rate constant. The standard deviation of the Gaussian noise was set at an adequate level (1% of the magnitude of each rate constant) to ensure proper coverage of search space while not sacrificing the optimization accuracy and search efficiency. (4) Steps 2 and 3 were repeated until all of the rate constants converged within acceptable ranges. The optimized rate constants were given as the averages of M particles. In Steps 1 and 2, the numerical calculation of the ODE was performed by using the backward differentiation formula (BDF) method. The program was parallelized by using OpenMP15 to improve computational performance. The CPU time needed to carry out the whole computation was 5 h on a workstation with two quad core Intel Xeon 2.8 GHz processors. Results Experimental Results. The pH oscillations measured under k0 ) 0.05 min-1 and [SO32-]0 ) 90 mM are illustrated in Figure 1. The oscillatory waveform exhibits typical relaxation oscillations characterized by the relatively slow pH change in both the low-pH and high-pH stages and rapid transitions between them. The minimum and maximum pH values are 3.7 and 6.7, respectively. We performed the model optimization based on the time series data in Figure 1. The experimental state diagram spanned by k0 and [SO32-]0 is shown in Figure 2. The system stays in the steady state with low pH (SSL) when [SO32-]0 is low, in which autocatalysis of H+ (D channel) is activated. On the other hand, if [SO32-]0 is high, the system stays in the steady state with high pH (SSH) where the autocatalysis is suppressed. In the conventional terminology, SSL is the thermodynamic branch, and SSH is the flow branch. The oscillatory state (OS) is observed in
Figure 2. Experimental state diagram in the k0-[SO23 ]0 plane. Fixed + 4parameters: [BrO3 ]0 ) 75 mM, [Fe(CN)6 ]0 ) 15 mM, [H ]0 ) 15 mM, and temperature ) 40.0 °C. The steady state with high/low pH (SSH/SSL) and the oscillatory state (OS) are indicated by closed up/ down triangles and open circles, respectively.
TABLE 3: List of the Rate Constants model
k3/M-1 s-1
5.97 × 10-2 original RKH model (35 °C)a optimized RKH 1.4 × 10-1 model (40 °C) extended RKH 9.8 × 10-2 model (40 °C) a
k4/M-1 s-1
k5/M s-1
k5′/M
18
1.5 × 10-5 2.5 × 10-4
13
9.6 × 10-6 1.0 × 10-3
22
2.9 × 10-5 8.0 × 10-4
Reference 12.
between the SSL and SSH. The obtained results are consistent with the preceding study by Edblom et al. (see Figure 5c in ref 11). The location of boundary lines is more clearly shown in the present state diagram by the higher-resolution experiment. Model Identification by the Particle Filter. In this subsection, we demonstrate the performance of the particle filter algorithm through the identification of a quantitatively accurate reaction model for the BSF system. We compare model calculations with experimental data for the three cases, (1) the RKH model with original rate constants (original RKH model), (2) the RKH model with optimized rate constants (optimized RKH model), and (3) the RKH model extended by the protonation equilibrium of SO42- with optimized rate constants (extended RKH model). (1) Original RKH Model. We first examined the properties of the original RKH model calculated with original rate constants (Tables 1 and 3). Figure 3 compares the calculated pH oscillations with the experimental data. Although the calculated pH waveform is qualitatively similar to the measured one, it has
Analysis of the BSF pH Oscillator
J. Phys. Chem. A, Vol. 114, No. 37, 2010 10093
Figure 3. The pH oscillations in the original RKH model compared with the experimental data. Rate constants are given in Table 3. Initial concentrations and k0 are the same as those in Figure 1.
Figure 6. State diagram obtained by the optimized RKH model. Fixed parameters and symbols are the same as those in Figure 2.
Figure 4. State diagram obtained by the original RKH model. Fixed parameters and symbols are the same as those in Figure 2.
Figure 7. The pH oscillations in the extended RKH model compared with the experimental data. Rate constants are given in Table 3. Initial concentrations and k0 are the same as those in Figure 1.
Figure 5. The pH oscillations in the optimized RKH model compared with the experimental data. Rate constants are given in Table 3. Initial concentrations and k0 are the same as those in Figure 1.
considerably larger amplitude and about a three times longer period. The state diagram obtained for the original RKH model is illustrated in Figure 4. The calculated boundary between OS and SSH shows a good agreement with the experiment (see Figure 2); however, the boundary between OS and SSL is displaced considerably to the low-[SO32-]0 and low-k0 side, indicating the insufficient stability of SSL in this model. (2) Optimized RKH Model. Next, we tried to optimize the RKH rate constants by using the particle filter without modifying the reaction scheme. The optimized rate constants are summarized in Table 3. Figure 5 shows the pH waveform calculated for the optimized RKH model. The oscillatory period becomes nearly the same as the measured one, with significantly reduced duration of the high-pH stage. This is mostly due to the substantial increase of k3. On the other hand, the calculated amplitude is still much larger than the observation. The state
diagram for the optimized RKH model is shown in Figure 6. As compared with the original RKH model, the stable region of SSL becomes slightly larger, whereas the SSH area is reduced, in accordance with the enhanced rate of reaction R3 and the retarded rate of reaction R5. At any rate, there exists a large discrepancy in the state diagram between the model and the experiment. (3) Extended RKH Model. Although the RKH model has succeeded in capturing the essential dynamical features of the BSF system, quantitative differences still remain between the model and the experiment even after optimizing the rate constants. Especially, a large discrepancy in the oscillatory minimum (of nearly 2 pH units) is not much improved. It is thus indicated that some additional processes are required for the model to obtain more quantitative agreement to the experiment. As the most probable assumption, we added the protonation equilibrium of SO42- as given in Table 2 (extended RKH model). The contribution of SO42- is considered to be important because SO42- is not only supplied externally in the present experiment but also generated internally through reactions R3 and R4.16 We applied the particle filter to the extended RKH model and reoptimized the rate constants. The obtained rate constants are summarized in Table 3. The calculated pH waveform illustrated in Figure 7 shows an excellent agreement to the experiment, including the level of oscillatory minimum. The state diagram for the extended RKH model is shown in Figure 8. The calculated boundary between OS and SSL is in good agreement with the experiment (see Figure 2). These improvements are naturally explained by the proton buffering effect induced by the additional equilibrium (reaction R6). The equilibrium reaction R6 suppresses excess increase of H+ during the low-pH stage in the oscillations. Furthermore, it also
10094
J. Phys. Chem. A, Vol. 114, No. 37, 2010
Figure 8. State diagram obtained by the extended RKH model. Fixed parameters and symbols are the same as those in Figure 2.
suppresses the rate of reaction R5 by limiting the concentration of free H+, which is considered to contribute to the stabilization of SSL. In Figure 8, we still recognize a small discrepancy with the experiment for the boundary location between OS and SSH; however, we may safely say that the extended RKH model with the optimized rate constants almost completely reproduces the experimental results. We also note that the optimized rate constants all lie within chemically acceptable ranges,14 which strongly supports the validity of the model. It has been demonstrated from the model calculation that fitting a single kinetic curve may lead to erroneous conclusions and that simultaneous fitting of multiple kinetic curves gives more accurate results.17 In order to confirm the accuracy of the present single-curve fitting, we have performed a simultaneous curve fitting by using the two different oscillatory time traces -1 and [SO2measured at [SO23 ]0 ) 90 mM and k0 ) 0.05 min 3 ]0 ) 80 mM and k0 ) 0.1 min-1. The objective function was defined as the nonweighted sum of RSS calculated over the interval T ) 2300 s for each data set, where the pH values were normalized with reference to the minimum and maximum values in each data set prior to the calculation. The fitting results were in fairly good agreement with the measured pH oscillations, as in the case of single-curve fitting (see Supporting Information). The estimated rate constants matched with the results of singlecurve fitting to within 9%, indicating that, for the analysis of the BSF system, fitting several cycles of oscillation provides sufficient information to give a reliable set of rate constants. Discussion The definition of the objective function, which is essentially problem-specific, is crucial for the successful application of the particle filter. In this study, a rather unconventional definition was employed for the objective function, namely, the RSS of pH summed up over fairly long time interval (eq 1). In the conventional use of the particle filter in estimating discrete time models, the objective function is calculated simply from the data at a single time step (usually formulated as a quadratic function), based on which a new particle set at the next time step is reproduced. By repeating the reproduction cycle, unsuitable particles are discarded, and new potential candidates are introduced around the suitable particles. In the present case of the BSF system, however, the instantaneous objective function such as the squared error of pH at a given time did not work, and no convergence was attained for the particle distribution. The problem is concerned with the drastic switching of the reaction channel inherent to the oscillatory chemical reactions. As with other chemical oscillators, the BSF system shows a sharp switching of the dominant reaction path due to the strong
Sato et al.
Figure 9. (a) Trends of the magnitude of pH sensitivities SipH (see eq 3 in the main text) during the oscillation cycle. (b) Reference pH oscillation. Calculation was performed by using the extended RKH model with the optimized rate constants (Table 3). Initial concentrations and k0 are the same as those in Figure 1. The sign of SipH is minus for i ) 3, 4, and 5′ and plus for i ) 5 throughout the oscillation.
nonlinearity induced by the autocatalysis (D channel). Accordingly, the contribution of each component reaction to the pH movement significantly varies as a function of the oscillatory phase, which strongly affects the evaluation of the objective function. The nonuniformity of pH response against the rate constant perturbation can be estimated by using the sensitivity analysis.18 We define the pH sensitivity as
SipH(t) )
∂pH(t) ∂ ln ki
for i ) 3, 4, 5, and 5′
(3)
which represents the magnitude of pH response caused by the small variation of (log-normalized) ki. Figure 9a shows the trends of |SpH i | (i ) 3, 4, 5, and 5′) plotted for one oscillatory period. In the high-pH stage, |S3pH| is extremely high, showing the overwhelming dominance of reaction R3 (S channel). During the autocatalytic phase, |S4pH| becomes greater than |S3pH|, corresponding to the switching of the dominant path from the S to D channel. In the low-pH stage, the pH response is mostly governed by reactions R4 and R5, while the contribution of reaction R3 is kept considerably low. As can be seen from Figure 9a, the dynamic range of |SpH i | is as large as 3 orders of magnitude for all i, which seriously disturbs smooth convergence of the particle distribution. If the instantaneous objective function is used, for example, the distribution of k5 quickly converges to the optimal value in the low-pH stage where reaction R5 is dominant, but it tends to diffuse in the high-pH stage where reaction R5 barely contributes. In the high-pH stage, the pH sensitivity with respect to k5 is so low that the particles with unsuitable k5 would survive in the next reproduction with high probability. There exist two possible solutions to work around the problem. One is to define the noninstantaneous objective function as given in eq 1 by taking the evaluation interval (T) longer than one oscillatory period so that all of the involved reactions can contribute to the objective function. The other one is to control the standard deviation of the additional Gaussian noise with reference to the trends of pH sensitivity. In this study, we adopted the former solution because of its algorithm simplicity and possible versatility to other nonlinear systems, though it may be slightly less advantageous in the computational efficiency. Figure 10 shows the evolution of marginal particle distributions during the optimization of the extended RKH model. As can be seen, all of the rate constants show smooth
Analysis of the BSF pH Oscillator
J. Phys. Chem. A, Vol. 114, No. 37, 2010 10095
Figure 10. Marginal distribution of (a-d) k3, (e-h) k4, (i-l) k5, and (m-p) k5′ after (a, e, i, m) 0 (initial distribution), (b, f, j, n) 5, (c, g, k, o) 10, and (d, h, l, p) 20 (final) reproductions, observed for the optimization of the extended RKH model.
TABLE 4: Average pH Sensitivity over a Whole Oscillation Cyclea
a
i
|SipH|
3 4 5 5′
6.0 × 10-1 6.9 × 10-2 4.0 × 10-2 1.6 × 10-2
Calculated from the data in Figure 9.
convergence, and the optimal values are reached after 20 reproduction steps. The convergence speed of the rate constants is in the order of k3 > k4 ≈ k5 > k5′, which is in good agreement with the order of the pH sensitivity averaged over one oscillation cycle (Table 4). Conclusion In summary, we analyzed the pH oscillations measured for the BSF system by using the particle filter to obtain a quantitatively accurate reaction model. The original RKH model was extended by adding the protonation equilibrium of SO42-. The extended RKH model with the optimized rate constants almost completely reproduced the measured pH waveform as well as the state diagram. Chemical oscillators such as the BSF system show drastic switching of the dominant reaction path, which makes it difficult to use the conventional objective function defined at a single time step. In this study, we proposed a noninstantaneous objective function with an evaluation interval longer than one oscillatory period, so that all of the relevant reaction steps can contribute to the objective function. Application of the present method to the analysis of related bromatesulfite-based pH oscillators19,20 is being planned in our laboratory in order to establish the generality of the present method and to promote a better quantitative understanding of these systems. The significance of the present study is that it has exemplified the availability of the particle filter algorithm to analyze real complex chemical systems. In particular, the reaction rate equations represented in the form of coupled differential
equations can be directly processed in the present method, which would be especially valuable for a chemist to perform prompt analysis of the reaction kinetics. It is expected that the present particle filter algorithm is applicable to the heuristic search of unidentified reaction mechanisms if combined with the appropriate reaction database system. The particle filter would also be effective for the biopathway analysis, where unidentified minor reactions may contribute as a significant source of nonGaussian noise. We are currently working on the development of the analysis software for more general use. Acknowledgment. N.S. and H.H.H. acknowledge N. Kumazawa and T. Washio for their suggestions, and support of this work. N.S. and H.H.H. thank T. Sugaya for his contributions. N.O. is obliged to Prof. Y. Cho (Tohoku University) for his kind encouragement. This work was supported by a Grantin-Aid for Scientific Research (Grant No.17540346) from the Ministry of Education, Culture, Sports, Science, and Technology of Japan. Supporting Information Available: Results of the twocurve fitting applied for the extended RKH model, including the table of rate constants, the calculated pH oscillations, and the state diagram. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Epstein, I. R.; Pojman, J. A. An Introduction to Nonlinear Chemical Dynamics; Oscillations, WaVes, Patterns, and Chaos; Oxford University Press: New York, 1998. (2) Kalman, R. E. J. Basic Eng. 1960, 82, 35–45. (3) Harvey, A. C. Forecasting, Structural Time Series Models and the Kalman Filter; Cambridge University Press: Cambridge, U.K., 1991. (4) Brown, R. G.; Hwang, P. Y. C. Introduction to Random Signals and Applied Kalman Filtering, 3rd ed.; Wiley: Hoboken, NJ, 1996. (5) Kitagawa, G. J. Comput. Graph. Stat. 1996, 5, 1–25. (6) Gordon, N. J.; Salmond, D. J.; Smith, A. F. M. IEE Proc. F: Commun., Radar, Signal Process. 1993, 140, 107–113. (7) Nakamura, K.; Yoshida, R.; Nagasaki, M.; Miyano, S.; Higuchi, T. Pac. Symp. Biocomput. 2009, 227–238.
10096
J. Phys. Chem. A, Vol. 114, No. 37, 2010
(8) Ng, B.; Pfeffer, A.; Dearden, R. Proceedings of the 19th International Joint Conference on Artificial Intelligence (IJCAI); Morgan Kaufmann Publishers: San Francisco, CA, 2005, pp 1360–1365. (9) Adachi, F.; Washio, T.; Motoda, H. Artificial Intelligence and Knowledge-Based Processing, IEICE Technical Report; 2004; Vol. 104, pp 23-28. (10) Epstein, I. R.; Kustin, K.; De Kepper, P.; Orbán, M. Sci. Am. 1983, 248, 112–123. (11) Edblom, E. C.; Luo, Y.; Orbán, M.; Kustin, K.; Epstein, I. R. J. Phys. Chem. 1989, 93, 2722–2727. (12) Rábai, G.; Kaminaga, A.; Hanazaki, I. J. Phys. Chem. 1996, 100, 16441–16442. (13) Hanazaki, I.; Rábai, G. J. Chem. Phys. 1996, 105, 9912–9920.
Sato et al. (14) Szirovicza, L.; Boga, E. Int. J. Chem. Kinet. 1998, 30, 869–874. (15) Openmp.org Homepage. http://www.openmp.org/ (2009). (16) Rábai, G.; Hanazaki, I. J. Phys. Chem. 1996, 100, 10615–10619. (17) Kormányos, B.; Horváth, A. K.; Peintler, G.; Nagypál, I. J. Phys. Chem. A 2007, 111, 8104–8109. (18) Turányi, T. J. Math. Chem. 1990, 5, 203–248. (19) Okazaki, N.; Rábai, G.; Hanazaki, I. J. Phys. Chem. A 1999, 103, 10915–10920. (20) Szántó, T. G.; Rábai, G. J. Phys. Chem. A 2005, 109, 5398– 5402.
JP106700F
