Revised dynamics of the Belousov-Zhabotinsky reaction model
aa r X i v : . [ n li n . C D ] M a y REVISED DYNAMICS OF THE BELOUSOV-ZHABOTINSKY REACTIONMODEL
JUDITA NAGYOV ´A, BRANISLAV JANS´IK, AND MAREK LAMPART
Abstract.
The main aim of this paper is to detect dynamical properties of the Gy¨orgyi-Fieldmodel of the Belousov-Zhabotinsky chemical reaction. The corresponding three-variable modelgiven as a set of nonlinear ordinary differential equations depends on one parameter, the flow rate.As certain values of this parameter can give rise to chaos, the analysis was performed in order toidentify different dynamics regimes. Dynamical properties were qualified and quantified using clas-sical and also new techniques. Namely, phase portraits, bifurcation diagrams, the Fourier spectraanalysis, the 0-1 test for chaos, and approximate entropy. The correlation between approximateentropy and the 0-1 test for chaos was observed and described in detail. Moreover, the three-stagesystem of nested subintervals of flow rates, for which in every level the 0-1 test for chaos andapproximate entropy was computed, is showing the same pattern. The study leads to an openproblem whether the set of flow rate parameters has Cantor like structure. Introduction
The Belousov-Zhabotinsky chemical reaction (BZ reaction), originally discovered in the 1950s byBoris P. Belousov [30], is an example of oscillating chemical reaction which can be maintained far fromequilibrium by an internal source of energy [2] resulting in a nonlinear chemical oscillator exhibitingdifferent dynamical regimes. Later on, the chemical mechanism of the reaction was described in [4],what is commonly called the FKN mechanism.There are many mathematical models representing different aspects of the BZ reaction. Forexample, the Brusselator, Oregonator and Gy¨orgyi-Field are three mathematical models for a typeof autocatalytic reaction – like the BZ reaction.The Oregonator model is the result of a quantitative kinetic analysis of oscillations in the BZreaction by Field and Noyes in 1974 [5] and it is a simplified version of the mechanism developed byField, K¨or¨os and Noyes (FKN mechanism) [30].The Brusselator model, a theoretical model for a type of autocatalytic reaction, was proposed byI. Prigogine and his collaborators [6].Finally, the Gy¨orgyi-Field model (GF model), describes the reaction taking place in a continuous-flow stirred-tank reactor (CSTR) [10], by a relatively simple mathematical model (see also [11] and[9]). This model, for a specific choice of parameters, exhibits chaos (see e.g. [3] and referencestherein or main results of this paper), on contrary to the Oregonator which has no chaotic solutions[31] describing the oscillatory behaviour and pattern formation in the BZ reaction. The GF modelwill be taken into consideration for further research in this paper.In recent decades, the BZ reaction has been extensively studied by physical chemists on its kineticbehaviour [3, 23, 29] and by mathematicians on the dynamics and patterns of the solutions of theassociated mathematical model [26, 20, 24, 31].More specifically, the research from the theory of dynamical system point of view was done. Thetransition from steady state to quasi-periodic and bursting oscillations, and further on to regularrelaxation oscillation via a complicated sequence of alternating periodic and chaotic regimes were done by simulations in [19]. The results of computer experiments on information processing in ahexagonal array of vesicles filled with BZ solution in a sub-excitable mode were introduced by [1].The discretized version of BZ reaction models was also researched. E.g. in [15], the dynamics ofthe local map is discussed, the set of trajectories that escape to infinity as well as analyze the set ofbounded trajectories – the Julia set of the system. The evidence of chaos was also demonstrated inan experimental way by dozens works e.g. [25, 32, 13].Despite a large number of results in the given area, it is possible to apply new methods to the givenBZ reaction model and to obtain new very interesting results that better characterize the trajectorybehaviour depending on the choice of state parameters showing properties of the parameter space.This work focuses on the characterization of dynamical properties of the GF model [10] dependingon the flow rate, denoted by k f . The qualitative and quantitative characterization of the dynamicsregimes is mainly done using the 0-1 test for chaos and approximate entropy. Recall that these toolswere applied in [16] to the two-dimensional coupled map lattice model of the Lagrangian type, whichis a discrete version of the BZ reaction. These tools were applied to the huge simulation data thatwas performed on the Salomon supercomputer at IT4Innovations National Supercomputing Centerlocated in Ostrava, Czech Republic.The paper is organized as follows: in Section 2 the model is introduced, followed by its mathemat-ical model in Section 3. The main results obtained by phase portraits, the Fourier spectra analysis,the approximate entropy, and the 0-1 for chaos, are contained in Section 4. Finally, the outcomesare summarized in Section 5.2. The Gy¨orgyi-Field reaction model
The GF model of the BZ reaction develops a description of the reaction in terms of a set ofdifferential equations containing only three variables. In common with experiments, the GF modelshows both regular, intermittent and chaotic behavior. While remaining close to a real chemicalsystem, it is sufficiently simple to allow detailed mathematical analysis [10]. The mechanism of thereaction is defined by the set of the following equations (1): Y + X + H → V (1a) Y + A + 2 H → V + H (1b) 2 X → V (1c) 0 . X + A + H → X + Z (1d) X + Z → . X (1e) V + Z → Y (1f) Z + M → (1g)where the corresponding chemical components are: Y = Br - , X = HBrO , Z = Ce , V =BrCH(COOH) or BrMA, A = BrO , H = H + , M = CH (COOH) . The concentrations of themain reactants A , H , M , and the total concentration of cerium ions C are summarized in Table 2. EVISED DYNAMICS OF THE BELOUSOV-ZHABOTINSKY REACTION MODEL 3
Table 1.
Rates and rate constants of the GF model chemical scheme.Reaction equation Rate r i Rate constant k i (1a) r = k HY X k = 4 . × M − s − (1b) r = k AH Y k = 2 . M − s − (1c) r = k X k = 3000 M − s − (1d) r = k A . H . ( C − Z ) X . k = 55 . M − . s − (1e) r = k XZ k = 7000 M − s − (1f) r = αk ZV k = 0 . M − s − (1g) r = βk M Z k = 0 . M − s − Mathematical model
A three-variable mathematical model of the BZ reaction, presented by Gy¨orgyi and Field in[10], describes the reaction taking place in a CSTR. The corresponding set of nonlinear ordinarydifferential equations contains only three variables, while still being able to accurately reproduce thebehavior of the BZ reaction observed experimentally [10] and it is based on a four-variable chemicalmechanism (1), see [10].The mathematical model, in its dimensionless form, consists of a set of scaled differential equa-tions: dxdτ = T ( − k HY x ˜ y + k AH Y /X ˜ y − k X x + 0 . k A . H . X − . ( C − Z z ) x . − (2a) − . k Z xz − k f x ) dzdτ = T (cid:0) k A . H . X . ( C/Z − z ) x . − k X xz − αk V zv − βk M z − k f z (cid:1) (2b) dvdτ = T (cid:0) k HX Y /V x ˜ y + k AH Y /V ˜ y + k X /V x − αk Z zv − k f v (cid:1) (2c)where τ = t/T , x = X/X , z = Z/Z , v = V /V , and ˜ y = (cid:0) αk Z V zv/ (cid:0) k HX x + k AH + k f (cid:1)(cid:1) /Y while t corresponding to time, X to HBrO , Y to Br - , Z to Ce , and V to BrMA. The rate constantsand parameters used in the following computations are given in the Table 1 and 2, respectively.The behavior of this system depends on the inverse residence time of the CSTR, the flow rate,noted k f [ s − ]. As certain values of this parameter can give rise to chaos, the following analysis wasperformed in order to identify different dynamics. Table 2.
Parameters of the investigated system (2).List of parameters A = 0 . M = 0 . H = 0 . C = 0 . α = 666 . β = 0 . J. NAGYOV´A, B. JANS´IK, AND M. LAMPART Main results
The system of differential equations 2 was solved numerically in Matlab [18] using the ode
45 solver.The simulations were done depending on free parameter k f ranging from 3 × − to 5 × − with10 − step. Each simulation was performed for the final time τ = 100 with a time step 10 − . Toavoid the systems distortions, only the last 20% of simulations were used for further computations.In all cases, initial conditions were set x = z = v = 1 . As a main results, phase diagrams, amplitude frequency spectrums (FFT), and Poincar´e sectionswere done for relevant choices of the parameter k f . To illustrate changes in dynamical behavior,bifurcation diagrams underlined by the approximate entropy and the 0-1 test for chaos with suitablemagnifications to the parameter k f were plotted.Consequently, as a goal of this paper, bifurcation diagrams, the approximate entropy, and the 0-1test for chaos were computed for nested set of parameters k f . The 0-1 test for chaos splits the valuesof the parameter for which regular (periodic or quasi-periodic) and irregular (chaotic) movementsappear, while the output of approximate entropy detects increasing complexity of the investigatedsystem 2.4.1. Phase diagrams, the Fourier spectra, and bifurcation diagrams.
Periodic as well aschaotic dynamics was identified in the studied model (2). For example, in Fig. 1 and Fig. 2, regularmovement is shown by a trivial loop (Fig. 1) for k f = 3 × − , and a non-trivial loop (Fig. 2) for k f = 3 . × − . Figure 3 gives an example of a chaotic trajectory, that is for k f = 3 . × − .The Fourier spectra were computed by Fast Fourier transform for k f = 3 × − , k f = 3 . × − ,and k f = 3 . × − , shown in Figs. 1, 2, and 3 respectively. Regular behavior is observable for thefirst two and chaos in the last case.In the case of regular movement, in Fig. 1 and Fig. 2, the Fourier spectra is formed by a numberof harmonic frequencies, hence the frequency of the periodic trajectory is computable. Periodicmotions of trajectory is also visible in Poincar´e sections.In the case of chaotic movement, Fig. 3, the Fourier spectra is formed by a number of harmoniccomponents having the basic, super-harmonic, sub-harmonic, and combination frequencies on whichfurther motions with frequencies forming the sided bands of the dominant frequencies are superposed.Their mutual ratio indicates the irregularity of the motion. The character of this chaotic motion isunderlined by the Poincar´e section.Next, the bifurcation diagram (constructed as a projection of a local maxima) of the model (2)was plotted for each variable x , z and v with respect to the free parameter k f ∈ (3 × − , × − ) Figure 1.
Phase diagram with Poincar´e section by plane v = 0 . x for k f = 3 × − . EVISED DYNAMICS OF THE BELOUSOV-ZHABOTINSKY REACTION MODEL 5
Figure 2.
Phase diagram with Poincar´e section by plane v = 0 . x for k f = 3 . × − . Figure 3.
Phase diagram with Poincar´e section by plane v = 0 .
847 and Fouriertransform of variable x for k f = 3 . × − .in Figure 4. In this bifurcation diagram, so-called ”period doubling” and ”windows” effects are alsovisible. Periodic movement can be identified in the range of the parameter, e.g., k f ∈ (3 × − , . × − ) and k f ∈ (3 . × − , × − ). The interval in between these values is interrupted by chaoticcases around k f = 3 .
25, and some k f ∈ (3 . × − , . × − ) and k f ∈ (3 . × − , . × − ). Figure 4.
Bifurcation diagram of the variable x (left), z (center), v (right) k f ∈ [3 × − , × − ].4.2. Approximate entropy.
The approximate entropy is a technique used to quantify the amountof regularity and unpredictability of fluctuations in time-series. The main advantages are that itcan be computed on short time series and it allows to compare the differences in complexity ofthe same system with different parameters settings, see, e.g., [21] to be detected. More complexnotion of entropy type can be found in, e.g., [14]. To compute the approximate entropy, two
J. NAGYOV´A, B. JANS´IK, AND M. LAMPART parameters must be set: embedding dimension m and neighborhood threshold r . Let s ( t ) ∈ R for t = { , , . . . , N } be a time series with N observations. Then embedded vector S ( t ) at time t , isdefined as S ( t ) = [ s ( t ) , s ( t + 1) , s ( t + 2) , . . . , s ( t + ( m − t is the observed time and m isthe embedding dimension. The maximum distance of embedded vectors is computed as follows:(3) D ( i, j ) = d ( S ( i ) , S ( j )) = max k =0 , ,...,m − | s ( i + k ) − s ( j + k ) | , for i, j = { , , . . . , N − ( m − } .Compute the thresholded version of the distance with threshold given by r :(4) d r ( i, j ) = (cid:26) , D ( i, j ) < r , otherwise, for i, j ∈ { , , . . . , N − ( m − } . Compute C mi ( r ) as a ratio between points in the neighborhood of i and the number of embeddedvectors.(5) C mi ( r ) = P N − ( m − j =1 d r ( i, j ) N − ( m − . Then compute the average of logarithm of all the C mi ( r )(6) Φ m ( r ) = 1 N − ( m − N − ( m − X i =1 ln C mi ( r ) . Finally, approximate entropy for the finite time series with N data points is computed as(7) ApEn ( m, r, N ) = Φ m ( r ) − Φ m +1 ( r ) . For robust estimation, it was suggested by Pincus [21] the time series is containing at least 10 observations.The approximate entropy was calculated using the TSEntropies package [28] for R [27]. Thecomputations were made for the input vector s given in a normalized form of all state variables: s ( t ) = p x ( t ) + z ( t ) + v ( t ) ,k f ∈ (3 × − , × − ) and r = 0 .
1. The results of approximate entropy for all values of theparameter k f are in Figs. 7–9.4.3. The 0-1 test for chaos, invented by Gottwald and Melbourne [7], is one ofthe methods for distinguishing between regular and chaotic dynamics of a deterministic system. Incontrast to the other approaches, the nature of the system is irrelevant, thus the test can be applieddirectly onto experimental data, ordinary differential equations, or partial differential equations.The results return values close to either 0 or 1, with 0 corresponding to regular dynamics and 1 tochaotic dynamics. With its easy implementation, evaluation, and wide range of application, usingthis tool for detecting chaos is becoming more popular in different fields.The 0-1 test for chaos can be computed by the following algorithm [7].Given the observation φ ( j ) for j = 1 , , ..., N and a suitable choice of c ∈ (0 , π ), the followingtranslation variables are computed: p c ( n ) = n X j =1 φ ( j ) cos( jc ) ,q c ( n ) = n X j =1 φ ( j ) sin( jc ) EVISED DYNAMICS OF THE BELOUSOV-ZHABOTINSKY REACTION MODEL 7
Figure 5.
The plot of p versus q for c = 1 . k f = 3 × − showingregular dynamics (left) and chaotic dynamics for k f = 3 . × − (right).for n = 1 , , ..., N . The dynamics of the translation components p c and q c is shown on the p c versus q c plot. A bounded trajectory is in Fig. 5 (left) corresponding to regular movement, for k f = 3 × − .An unbounded trajectory is in Fig. 5 (right) related to the chaotic case, for k f = 3 . × − .The idea for the 0-1 test, first described in [7], is that the boundedness or unboundedness of thetrajectory { ( p j , q j ) j ∈ [1 ,N ] } can be studied through the asymptotic growth rate of its time-averagedmean square displacement (MSD), which is defined as M ( n ) = lim N →∞ N N X j =1 d ( j, n ) where d ( j, n ) = q ( p j + n − p j ) + ( q j + n − q j ) is the time lapse of the duration n ( n ≪ N ) starting from the position at time j . As it is shownin [7, 8], it is important to use values of n small enough compared to N , noted n cut , ( n ≤ n cut ). Asubset of time lags n cut ∈ [1 , N/
10] is advised for the computation of each K c .For bounded trajectories and regular dynamics, M ( n ) is a bounded function in time, whereasunbounded trajectories, meaning chaotic dynamics, are described by M ( n ) growing linearly withtime. Thus the asymptotic growth rate of the MSD must be calculated, which correlates with theunboundedness of the trajectory.As proposed in [7], the modified MSD is calculated as D ( n ) = M ( n ) − E ( φ ) − cos( nc )1 − cos c The output of the 0-1 test for chaos is computed by the correlation method as K c = corr( ξ, ∆) ∈ [ − , ξ = (1 , , ..., n cut ) and ∆ = ( D c (1) , D c (2) , ..., D c ( n cut )).The final result of the test is K = median( K c ) . J. NAGYOV´A, B. JANS´IK, AND M. LAMPART
Figure 6.
The plot of K c based on c for k f = 3 × − showing regular dynamics(left) and for k f = 3 . × − showing chaotic dynamics (right). Figure 7.
The result of approximate entropy (in blue) and the result of the 0-1test for chaos (in red) for k f ∈ (3 × − , × − ). The magnification in the blackrectangle is shown in Fig 8. The bifurcation diagram for variable x is shown in thebackground.The position of the studied system (2) at any moment of time is determined by displacements x , z , and v , which are used for defining vector φ : φ ( j ) = p x ( j ) + z ( j ) + v ( j ) . For these simulations, a free software environment R [27] was used including the
Chaos01 packagedeveloped by T. Martinoviˇc [17]. Comparison of values K c for periodic and chaotic case is shown inFig. 6, for k f = 3 × − and k f = 3 . × − , respectively.The results of the 0-1 test for chaos for all values of the parameter k f are in Figs. 7–9. EVISED DYNAMICS OF THE BELOUSOV-ZHABOTINSKY REACTION MODEL 9
Figure 8.
The result of approximate entropy (in blue) and the result of the 0-1test for chaos (in red) for k f ∈ (3 . × − , . × − ). The magnification in theblack rectangle is shown in Fig 9. The bifurcation diagram for variable x is shownin the background. Figure 9.
The result of approximate entropy (in blue) and the result of the 0-1test for chaos (in red) for k f ∈ (3 . × − , . × − ). The bifurcation diagramfor variable x is shown in the background.5. Conclusions
In this paper, the GF model (2) associated with the BZ chemical reaction was solved usingadaptive six-stage, fifth-order, Runge-Kutta method implemented as ode
45 solver in Matlab. Toeliminate the stiffness problem the model (2) was also simulated by ode s solver in Matlab [18],outputs were identical.The simulations were used to plot 3D phase portraits, bifurcation diagrams, the approximateentropy, and the 0-1 test for chaos. The mining process of dynamical properties were performed inthe free software R [27] using packages T SEntopies [28] and
Chaos
01 [17] depending on flow rateparameter k f .It is evident from main results in Figs. 7–9 that both tests clearly detect regular and irregularpatterns for given k f . Our results show that the method of approximate entropy returns qualification constant whichdescribes complexity in the system invariantly with respect to the origin. On the other hand, the0-1 test as qualification tool returns zero for regular (periodic or quasi-periodic) movement and onefor irregular (chaos) characteristics.Moreover, if the output of the 0-1 test is not close to zero or one, then the examined test case hasnot yet reached attractor or reached intermittent state, see e.g. [12] or [22] and references therein.Further, we observe a correlation between approximate entropy and the 0-1 test for chaos. Ingeneral, the increasing values of the 0-1 test for chaos are coupled to increasing approximate entropyand vice versa.We notice isolated low values of the 0-1 test for chaos accompanied by comparatively low valuesof approximate entropy well within of chaotic region characterized by high 0-1 test chaos valuesand approximate entropy. To investigate and zoom in, we have constructed three-stage system ofnested subintervals of flow rates k f , see Figs. 7–9, for which in every level the 0-1 test for chaos andapproximate entropy was computed. At every level we observe the same pattern.That naturally yields suggestion of a fractal structure in the set of k f : Open Problem 1.
Is there a totally disconnected (Cantor) set of flow rates k f in [3 × − , × − ] such that for each such parameter the GF model (2) is showing chaos? Open Problem 2.
Is there a totally disconnected (Cantor) set of flow rates k f in [3 × − , × − ] such that for each such parameter the GF model (2) is showing regular movement? Acknowledgement
This work was supported by The Ministry of Education, Youth and Sports from the NationalProgramme of Sustainability (NPU II) project “IT4Innovations excellence in science – LQ1602“;by The Ministry of Education, Youth and Sports from the Large Infrastructures for Research, Ex-perimental Development and Innovations project “IT4Innovations National Supercomputing Center– LM2015070“; by SGC grant No. SP2019/125 ”Qualification and quantification tools applicationto dynamical systems”, VˇSB - Technical University of Ostrava, Czech Republic, Grant of SGS No.SP2019/84, VˇSB - Technical University of Ostrava, Czech Republic.
References [1] A. Adamatzky, J. Holley, L. Bull, B. De Lacy Costello, Chaos Solitons Fractals 44, 779 (2011)[2] R. D’Ambrosio, M. Moccaldi, B. Paternoster, F. Rossi, J. Math. Chem. 56, 2876 (2018)[3] I.R. Epstein, J.A. Pojman, An Introduction to Nonlinear Chemical Dynamics: Oscillations, Waves, Patterns, andChaos (Oxford University Press, 1998)[4] R.J. Field, E. K¨or¨os, R.M. Noyes, J. Am. Chem. Soc. 94, 8649 (1972)[5] R.J. Field, R.M. Noyes, J. Chem. Phys. 60, 1877 (1974)[6] P. Glansdorff, I. Prigogine, Thermodynamic Theory of Structure, Stability and Fluctuations (Wiley, 1971)[7] G.A. Gottwald, I. Melbourne, SIAM J. Appl. Dyn. Syst. 8, 129 (2009)[8] G.A. Gottwald, I. Melbourne, Nonlinearity 22, 1367 (2009)[9] L. Gy¨orgyi, R.J. Field, J. Phys. Chem. 95, 6594 (1991)[10] L. Gy¨orgyi, R.J. Field, Nature 355, 808 (1992)[11] L. Gy¨orgyi, S. L. Rempe, R. J. Field, J. Phys. Chem. 95, 3159 (1991)[12] J.F. Heagy, N. Platt, S.M. Hammel, Phys. Rev. E 49, 1140 (1994)[13] J.L. Hudson, M. Hart, D. Marinko, J. Chem. Phys. 71, 1601 (1979)[14] W. Jun, T. Lingyu, Z. Xianyong, L. Yuyan, Appl. Math. Nonlinear Sci. 2, 329 (2017)[15] H. Kang, Y. Pesin, Milan J. Math. 73, 1 (2005)[16] M. Lampart, T. Martinoviˇc, J. Math. Chem. 57, 1670 (2019)[17] T. Martinoviˇc, Chaos01: 0-1 Test for Chaos. (2016), https://CRAN.R-project.org/package=Chaos01. Accessed1 October 2018[18] MATLAB 2015b (The MathWorks, Inc., 2015)
EVISED DYNAMICS OF THE BELOUSOV-ZHABOTINSKY REACTION MODEL 11 [19] O.V. Noskov, A.D. Karavaev, V.P. Kazakov, S.I. Spivak, Mendeleev Commun. 4, 82 (1994)[20] C.V. Pao, J. Part. Diff. Eq. 1, 61 (1988)[21] S.M. Pincus, Proc. Natl. Acad. Sci. U.S.A. 88, 2297 (1991)[22] N. Platt, E.A. Spiegel, C. Tresser, Phys. Rev. Lett. 70, 279 (1993)[23] I. Prigogine, R. Lefever, J. Chem. Phys. 48, 1695 (1968)[24] W.H. Ruan, C.V. Pao, J. Math. Anal. Appl. 169, 157 (1992)[25] R.A. Schmitz, K.R. Graziani, J.L. Hudson, J. Chem. Phys. 67, 3040 (1977)[26] F. Schneider, Angew. Chem. 98, 941 (1986)[27] R: A Language and Environment for Statistical Computing (R Core Team, R Foundation for Statistical Com-puting, 2018)[28] J. Tomˇcala, TSEntropies: Time Series Entropies. (2018), https://CRAN.R-project.org/package=TSEntropies.Accessed 22 March 2019[29] J.J. Tyson, J. Phys. Chem. 86, 3006 (1982)[30] J.J. Tyson, in Lecture Notes in Biomathematics (Springer, 1994), p. 569[31] J.J. Tyson, P.C. Fife, J. Chem. Phys. 73, 2224 (1980)[32] H. Yamazaki, Y. Oono, K. Hirakawa, J. Phys. Soc. Jpn. 46, 721 (1979)
IT4Innovations, VˇSB - Technical University of Ostrava, 17. listopadu 15/2172, 708 33 Ostrava, CzechRepublic, Department of Applied Mathematics, VˇSB - Technical University of Ostrava, 17. listopadu15/2172, 708 33 Ostrava, Czech Republic
E-mail address : [email protected] IT4Innovations, VˇSB - Technical University of Ostrava, 17. listopadu 15/2172, 708 33 Ostrava, CzechRepublic
E-mail address : [email protected] IT4Innovations, VˇSB - Technical University of Ostrava, 17. listopadu 15/2172, 708 33 Ostrava, CzechRepublic, Department of Applied Mathematics, VˇSB - Technical University of Ostrava, 17. listopadu15/2172, 708 33 Ostrava, Czech Republic
E-mail address ::