Experimentally observable transitions between dynamical states in complex reaction systems
Introduction
Complex reaction systems are very common in nature, as well as, in human made environments (Nicolis & Prigogine, 1990) Atmosphere chemistry (Clapp & Jenkin, 2001; Finlayson-Pitts & Pitts, 1997; Prinn, 2003) and realistic industrial catalytic processes (Berezin, Denisov, & Emanuel, 1969; Coppens & Froment, 1996; Quiceno, Perez-Ramirez, Warnatz, & Deutschmann, 2006) are examples of very complex systems. Also, biochemical reaction networks with hundreds of different reaction species and elementary reactions are subject of recent investigation (Cox et al., 2006; Palsson, Price, & Papin, 2003) Without a doubt, in complex reaction systems one can found fascinating kinetic phenomena, like multi-stability, excitability, oscillations, traveling waves and chaos (Anić et al., 1998, Field and Burger, 1985; Gray & Scott, 1990; Kolar-Anić, Stanisavljev, Krnajski Belovljev, Peeters, & Anić, 1990; Schmitz, 1990; Schmitz, Kolar-Anić, Anić, & Čupić, 2000; Schmitz, Kolar-Anić, Anić, Grozdić, & Vukojević, 2006; Scott, 1991). Moreover, variety of different shapes in the signal time sequences may occur in almost any chemical reaction, during the transition periods. Any reaction system starting from some arbitrary initial conditions and traveling to final steady or equilibrium state has to pass through the distant regions of phase space while number of different chemical species is simultaneously transformed through the reaction network (Higgins, 1965, Schmitz, 1990).
On the other hand, it is rarely necessary to know all details of the reaction mechanism of most processes proceeding under well controlled and only slightly variable reaction conditions. In most cases, reaction dynamics is reduced to lower dimensional attractors, simplifying the underlying mathematical models. According to (Higgins, 1965), who did early systematic study of the oscillatory reactions (Higgins, 1967), especially in biochemical systems, the reaction dynamics, in general, can be reduced to transitions between different dynamical states. In order to illustrate the changes in the reaction system regimes Higgins developed a concept of dynamical state level diagrams, based on the general description of the dynamical systems with different time scales (Fig. 1). Motion of the dynamical system in phase space, was depicted in Higgins's work through the transition between distinct states – special states in the state space, where the motion of the dynamical system obeys a reduced set of equations – dynamical states. Moreover, system moves from one to another dynamical state within a characteristic time, more specifically half-life period of the dynamical state.
Higgins (1965) used terms of new species, defined as equilibrated fast species. In the Higgins's paper (1965) this new hypothetical species corresponds to some specific dynamical state of the system. For example, in the case of two consecutive reversible reactions illustrated in Fig. 1, the symbol A(eq)B stands for the new hypothetical species, corresponding to equilibrated A and B, and A(eq)B(eq)C corresponds to the equilibrium.
Starting from the arbitrary initial state , reaction system is subjected first to a very fast transition to new dynamical state and then also to another transition to its final equilibrium state A(eq)B(eq)C.
Higgins used dynamical state transitions to explain the control mechanism in cellular reactions with complex reaction kinetics, and further to develop the theory of oscillatory reactions (Higgins, 1967). He defined three kinds of species which differ in respect to the level of the corresponding dynamical state. Variables “well below” the level of the actual dynamical state (at frequency scale) can be approximated as constant and classified as structural control variables (SCV). Variables “well above” given level can be referred as hidden variables, since they will not generally be excited and they have no effect. These are the fast, dependent variables which do not enter reduced differential equations. Finally, between these extremes there are independent dynamical variables (Higgins, 1965, Higgins, 1967). Obviously, the choice of the dynamic variables is dependent on actual dynamical state and further on the time scale of interest.
Although the concept of dynamical state level diagrams is not widely accepted, reaction systems with multiple time scales are studied extensively. Schmitz defined criterion for quasi-steady states, that can also be useful for choice of dynamical variables and hence, for identification of actual dynamical states (Schmitz, 1983). It is based on the ratio between the values of net-flux and one-direction-flux through the particular species of the reaction network. In the equilibrium state net fluxes are identical to zero for all species, since the rates of formation are equal to the rates of the destruction:where r+ is the rate of formation, and r− is the rate of destruction for some species.
However, in non-equilibrium systems there is always non-zero net reaction in which some reactants are spent and some products are formed. Therefore, the rate of concentration changes of these external species is also non-zero. There is usually other kind of species which are formed and destroyed by the reaction network processes, and they will be called internal ones. The rate of change of the internal species concentration is then:where the summation is over all the reactions in which the specified species are formed or destroyed, respectively. If internal species is formed and destroyed with equal rates, steady state approximation is applicable:
If the reaction rates are dependent on slowly changing concentrations of the external species (reactants and products), then the steady state concentrations of the internal species also depend on time. Therefore, the rates of formation or destruction of the internal species will not be exactly equal and we have quasi-steady-state approximation:The choice of internal species in this case depends on the time scale of interest. According to Schmitz criterion (Schmitz, 1983), we can say that some species is in quasi-steady-state if the rate of its concentration change is much lower than both, the sum of the rates of its formation and the sum of the rates of its destruction:
We can use Schmitz criterion (5) to define which species are the equilibrated ones, according to Higgins terminology, and hence to identify dynamical state of the system.
Also, other forms of transitions between different dynamical states were studied in the multiple time scale systems (Gray & Scott, 1990). It seems that the transition between dynamical states or quasi-steady states in chemical reaction systems is connected with changes in dominating reaction pathways. Moreover, different reaction pathways can sometimes be ascribed to different kinetic behavior of the reaction systems at different time scales. Bray already noticed in the first description of a homogeneous oscillatory reaction (Bray, 1921), that during the oscillations, sequence of periods can be observed, where the domination is exchanged between the reduction of iodate and the oxidation of iodine by hydrogen peroxide. This oscillatory reaction, now known as the Bray–Liebhafsky (BL) reaction, is hydrogen peroxide decomposition to water and oxygen with iodate ions as a catalyst, present in acid medium (Anić, Mitić, & Kolar-Anić, 1985; Bray, 1921; Bray & Liebhafsky, 1931; Matsuzaki, Alezander, & Liebhafsky, 1970; Woodson & Liebhafsky, 1969). In their papers (Kolar-Anić, Čupić, Anić, & Schmitz, 1997; Schmitz, 1998) the authors described BL oscillatory reaction dynamics as a dynamical system with multiple time scales. The progress of the reaction system in state space during the oscillations was described as a sequence of changes in dominating reaction pathways. Different reaction pathways in oscillatory reactions were also studied by Clark, Schmitz, Schreiber, Anić, Kolar-Anić, using the Stoichiometric Network Analysis (Clarke, 1980, Clarke, 1988, Schmitz, 1991; Schreiber, Hung, & Ross, 1996; Strasser, Stemwedel, & Ross, 1993). In further investigation of the BL system, activation energies were evaluated, corresponding to different reaction pathways and characterizing different segments of time series data (Anić, 1997, Anić et al., 1997; Ćirić, Anić, Čupić, & Kolar-Anić, 2000).
According to previous considerations it follows that dynamical systems with multiple time scales obey different kinetic laws in different dynamical states. This property was used in (Forsythe & Mavrovouniotis, 1997) to develop a special technique that enabled the authors to simplify the modeling of the reaction systems behavior passing in time from one to another region of the composition space. The inherent structure of the reaction systems was used in (Forsythe & Mavrovouniotis, 1997) to divide the composition space into regions. Within each region, the order-of-magnitude relationships between terms in rate equations are used systematically to reduce the order and coupling of the model. The full system was then described through a piecewise combination of these simpler, region-specific models.
It remains an open question at what level it is possible to identify dynamical state transitions from kinetic experiments. Piecewise descriptions of time series from kinetic experiments are a subject of different forms of trend analysis. In trend analysis, time intervals (episodes) are identified with identical combination of first and second derivative signs of the time concentration profiles (Table 1). This technique is found to be useful in numerous applications from chemical engineering to fast biochemical tests of blood samples taken from wounded soldiers on a battlefield (Charbonnier, Garcia-Beltan, Cadet, & Gentil, 2005; Dash, Maurya, Venkatasubramanian, & Rengaswamy, 2004; Melek, Lu, Kapps, & Fraser, 2005; Schaich, Becker, & King, 2001; Zaliapin, Gabrielov, & Keilis-Borok, 2004).
Table 1 would be more complete if two more combinations of derivative signs would be added: “0,+” and “0,−”, corresponding to kinetic curves with minimum and maximum values. However, it is usually more convenient to divide kinetic curve in separate trend episodes ending at these easily identifiable points.
It can be postulated that different trends are somehow connected with some sort of dynamical states. Due to possible importance of trend analysis for the detection of dynamical states and for the detection of transitions between them, a brief review of recent achievements in the area is given below.
Various approaches were used for trend detection in time series, from simple regression to fuzzy logic (Melek et al., 2005). Several problems arose in application of trend analysis, such as, numerical problems in treating noisy data, non-linearity detection, or model identification.
On-line extraction of semi-qualitative temporal episodes from any univariate time-series was developed in (Charbonnier et al., 2005). Only three primitives were used here to describe the episodes: Increasing, Decreasing and Steady. The method acts on noisy data, without prefiltering. It consists of several steps, including: (1) on-line segmentation of data into linear segments, (2) classification of the last segment into semi-quantitative episode, using three primitives, and (3) aggregation of the current episode to the previous one to form the trend. The method was designed to be useful engineering tool in control of the end product quality, depending on variations in raw materials as input (Charbonnier et al., 2005).
Another simple but effective approach was proposed (Zaliapin et al., 2004), based on piecewise linear approximation of time series. Computationally effective algorithm was constructed for the decomposition of a time series into a hierarchy of trends at different scales. This approach was particularly appropriate for self-affine analysis and obtaining of local fractal properties of the processes. It also allows detection of the non-linear correlations in non-linear, long-term trends (Zaliapin et al., 2004).
Automatic detection of trends was proposed (Dash et al., 2004) based on polynomial-fit based interval-halving technique, with wavelet based denoising. The least-order (among constant, first-order and quadratic) polynomial with fit-error statistically insignificant compared to noise (as dictated by F-test) is used to represent the segment. If the fit error is large even for the quadratic polynomial, then the length is halved and the process is repeated on the segments. Finally, a unique assignment of qualitative shape is made to each of the identified segments (Dash et al., 2004).
Finally, trend analysis was also the underlying method in a tool (TAM-C) developed for the automatic identification of mathematical models of chemical reaction systems (Schaich et al., 2001). The software imitates human ability to reason about system structure and system behavior at a qualitative level. Temporal shape of a continuous variable was described by a sequence of, so called, episodes, and this sequence is called history. Episode is a time interval in which signs of the first and the second time derivative of the variable are not changed. Each episode is described by the +, 0 or − sign combination for the two derivatives.
To overcome measurement noise which is usually present in real data, the obtained time series is smoothed by cubic splines here (Schaich et al., 2001). At the end, the authors sought for a model which gives the same qualitative history with transitions somewhere in the intervals determined by the experiments. Therefore, the so-called, interacting interval simulation was used here, for the interval identification in specific cases. Finally, the procedure was used to identify model parameters of the industrial process (Schaich et al., 2001).
In the present paper, identification of trends and their correspondence with dynamical state transitions are briefly discussed on several reaction network models. The influence of different threshold values on the trend analysis and dynamical state analysis results is examined. Both methods are applied here to the results of the numerical simulation, because, only in this case the actual dynamical states of the system are known. For this purpose the models of three reaction systems are used: simple model of consecutive, first order, equilibrium reactions, already discussed by Higgins (1965), partial oxidation of cyclohexane (Khar’kova, Arest-Yakubovich, & Lipes, 1989), and model M(1-8) of the oscillatory reaction Bray–Liebhafsky (Kolar-Anić, Mišljenović, Anić, & Nicolis, 1995a; Kolar-Anić, Vukelić, Mišljenović, & Anić, 1995b).
Section snippets
Numerical simulations
The numerical simulation was performed with the Gear algorithm (Gear, 1971), suitable for the integration of stiff differential equations.
Dynamical state analysis
Dynamical states were determined for every point of the time sequence, according to the slightly changed Schmitz criterion function:where ɛ is a threshold value. We found criterion (6) more efficient in computing then (5), since it contains only one inequality. It is also applicable in systems arbitrarily far from equilibrium, either
Model 1
Rate constants and initial concentrations in the Model 1 (equations in Section 2.4.1) simulation were chosen so that maximum occurred in concentration of intermediary B, after short time and dynamical system remains in quasi-steady state for a long time, as it is illustrated in Fig. 2, corresponding to the case when k1 + k−1 ≫ k2 + k−2 in Higgins's paper (1965). However, picture looks different on different time scales (Fig. 2).
The cross points with several threshold values used in dynamical state
Conclusions
The existence of the qualitative correspondence between dynamical state transitions and sequences of trend episodes is generally confirmed on all three analyzed reaction network models. Hence, it was clearly shown in the most complex case of the oscillatory reaction BL, Model 3, how transient dynamical states may be identified through the symbolic representation of the sequence of trend episodes in concentration time series. Application of the Stoichiometric Network Analysis on different
Acknowledgements
This work was supported by the Serbian Ministry of Science and Environmental Protection through the Project 142019B: The preparation, characterization and tests of the catalytic properties of the tailored materials.
References (60)
- et al.
Analysis of the relationship between ambient levels of O3, NO2 and NO as a function of NOx in the UK
Atmospheric Environment
(2001) - et al.
Trends extraction and analysis for complex system monitoring and decision support
Engineering Applications of Artificial Intelligence
(2005) - et al.
Fractal aspects in the catalytic reforming of naphtha
Chemical Engineering Science
(1996) - et al.
Development of a metabolic network design and optimization framework incorporating implementation constraints: A succinate production case study
Metabolic Engineering
(2006) Dynamics and control in cellular reactions
- et al.
The first maximum of the iodide concentrations in the Bray–Liebhafsky reaction
Computers Chemistry
(1990) - et al.
Development of network-based pathway definitions: the need to analyze real metabolic networks
Trends in Biotechnology
(2003) - et al.
Kinetic model of cyclohexane oxidation
Chemical Engineering Science
(2001) - et al.
Modeling the high-temperature catalytic partial oxidation of methane over platinum gauze: Detailed gas-phase and surface chemistries coupled with 3D flow field simulations
Applied Catalysis A: General
(2006) - et al.
Qualitative modelling for automatic identification of mathematical models of chemical reaction systems
Control Engineering Practice
(2001)