Revised dynamics of the Belousov-Zhabotinsky reaction model (original) (raw)

REVISED DYNAMICS OF THE BELOUSOV-ZHABOTINSKY REACTION MODEL

JUDITA NAGYOVÁ, BRANISLAV JANSÍK, AND MAREK LAMPART

Abstract

The main aim of this paper is to detect dynamical properties of the Györgyi-Field model of the Belousov-Zhabotinsky chemical reaction. The corresponding three-variable model given 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 to identify different dynamics regimes. Dynamical properties were qualified and quantified using classical and also new techniques. Namely, phase portraits, bifurcation diagrams, the Fourier spectra analysis, the 0−10-1 test for chaos, and approximate entropy. The correlation between approximate entropy and the 0−10-1 test for chaos was observed and described in detail. Moreover, the three-stage system of nested subintervals of flow rates, for which in every level the 0−10-1 test for chaos and approximate entropy was computed, is showing the same pattern. The study leads to an open problem whether the set of flow rate parameters has Cantor like structure.

1. InTRODUCTION

The Belousov-Zhabotinsky chemical reaction (BZ reaction), originally discovered in the 1950s by Boris P. Belousov [30], is an example of oscillating chemical reaction which can be maintained far from equilibrium by an internal source of energy [2] resulting in a nonlinear chemical oscillator exhibiting different 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. For example, the Brusselator, Oregonator and Györgyi-Field are three mathematical models for a type of autocatalytic reaction - like the BZ reaction.

The Oregonator model is the result of a quantitative kinetic analysis of oscillations in the BZ reaction by Field and Noyes in 1974 [5] and it is a simplified version of the mechanism developed by Field, Körös and Noyes (FKN mechanism) [30].

The Brusselator model, a theoretical model for a type of autocatalytic reaction, was proposed by I. Prigogine and his collaborators [6].

Finally, the Györgyi-Field model (GF model), describes the reaction taking place in a continuousflow 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 references therein 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 model will 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 kinetic behaviour [3,23,29][3,23,29] and by mathematicians on the dynamics and patterns of the solutions of the associated mathematical model [26,20,24,31][26,20,24,31].

More specifically, the research from the theory of dynamical system point of view was done. The transition from steady state to quasi-periodic and bursting oscillations, and further on to regular relaxation 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 a hexagonal 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 of the local map is discussed, the set of trajectories that escape to infinity as well as analyze the set of bounded trajectories - the Julia set of the system. The evidence of chaos was also demonstrated in an 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 given BZ reaction model and to obtain new very interesting results that better characterize the trajectory behaviour 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] depending on the flow rate, denoted by kfk_{f}. The qualitative and quantitative characterization of the dynamics regimes is mainly done using the 0−10-1 test for chaos and approximate entropy. Recall that these tools were applied in [16] to the two-dimensional coupled map lattice model of the Lagrangian type, which is a discrete version of the BZ reaction. These tools were applied to the huge simulation data that was performed on the Salomon supercomputer at IT4Innovations National Supercomputing Center located in Ostrava, Czech Republic.

The paper is organized as follows: in Section 2 the model is introduced, followed by its mathematical model in Section 3. The main results obtained by phase portraits, the Fourier spectra analysis, the approximate entropy, and the 0−10-1 for chaos, are contained in Section 4. Finally, the outcomes are summarized in Section 5.

2. The Györgyi-Field Reaction model

The GF model of the BZ reaction develops a description of the reaction in terms of a set of differential equations containing only three variables. In common with experiments, the GF model shows both regular, intermittent and chaotic behavior. While remaining close to a real chemical system, it is sufficiently simple to allow detailed mathematical analysis [10]. The mechanism of the reaction is defined by the set of the following equations (1):

Y+X+H→2VY+A+2H→V+H2X→V0.5X+A+H→X+ZX+Z→0.5XV+Z→YZ+M→\begin{aligned} Y+X+H & \rightarrow 2 V \\ Y+A+2 H & \rightarrow V+H \\ 2 X & \rightarrow V \\ 0.5 X+A+H & \rightarrow X+Z \\ X+Z & \rightarrow 0.5 X \\ V+Z & \rightarrow Y \\ Z+M & \rightarrow \end{aligned}

where the corresponding chemical components are: Y=Br−,X=HBrO2,Z=Ce4+,V=Y=\mathrm{Br}^{-}, X=\mathrm{HBrO}_{2}, Z=\mathrm{Ce}^{4+}, V= BrCH⁡(COOH)2\operatorname{BrCH}(\mathrm{COOH})_{2} or BrMA,A=BrO3−,H=H+,M=CH2(COOH)2\mathrm{BrMA}, A=\mathrm{BrO}_{3}{ }^{-}, H=\mathrm{H}^{+}, M=\mathrm{CH}_{2}(\mathrm{COOH})_{2}. The concentrations of the main reactants A,H,MA, H, M, and the total concentration of cerium ions CC are summarized in Table 2.

Table 1. Rates and rate constants of the GF model chemical scheme.

3. Mathematical model

A three-variable mathematical model of the BZ reaction, presented by Györgyi and Field in [10], describes the reaction taking place in a CSTR. The corresponding set of nonlinear ordinary differential equations contains only three variables, while still being able to accurately reproduce the behavior of the BZ reaction observed experimentally [10] and it is based on a four-variable chemical mechanism (1), see [10].

The mathematical model, in its dimensionless form, consists of a set of scaled differential equations:

dxdτ=T0(−k1HY0xy~+k2AH2Y0/X0y~−2k3X0x2+0.5k4A0.5H1.5X0−0.5(C−Z0z)x0.5−−0.5k5Z0xz−kfx)dzdτ=T0(k4A0.5H1.5X0.5(C/Z0−z)x0.5−k5X0xz−αk6V0zv−βk7Mz−kfz)dvdτ=T0(2k1HX0Y0/V0xy~+k2AH2Y0/V0y~+k3X02/V0x2−αk6Z0zv−kfv)\begin{aligned} \frac{d x}{d \tau}= & T_{0}\left(-k_{1} H Y_{0} x \tilde{y}+k_{2} A H^{2} Y_{0} / X_{0} \tilde{y}-2 k_{3} X_{0} x^{2}+0.5 k_{4} A^{0.5} H^{1.5} X_{0}^{-0.5}\left(C-Z_{0} z\right) x^{0.5}-\right. \\ & \left.-0.5 k_{5} Z_{0} x z-k_{f} x\right) \\ \frac{d z}{d \tau}= & T_{0}\left(k_{4} A^{0.5} H^{1.5} X^{0.5}\left(C / Z_{0}-z\right) x^{0.5}-k_{5} X_{0} x z-\alpha k_{6} V_{0} z v-\beta k_{7} M z-k_{f} z\right) \\ \frac{d v}{d \tau}= & T_{0}\left(2 k_{1} H X_{0} Y_{0} / V_{0} x \tilde{y}+k_{2} A H^{2} Y_{0} / V_{0} \tilde{y}+k_{3} X_{0}^{2} / V_{0} x^{2}-\alpha k_{6} Z_{0} z v-k_{f} v\right) \end{aligned}

where τ=t/T0,x=X/X0,z=Z/Z0,v=V/V0\tau=t / T_{0}, x=X / X_{0}, z=Z / Z_{0}, v=V / V_{0}, and y~=(αk6Z0V0zv/(k1HX0x+k2AH2+kf))/Y0\tilde{y}=\left(\alpha k_{6} Z_{0} V_{0} z v /\left(k_{1} H X_{0} x+k_{2} A H^{2}+k_{f}\right)\right) / Y_{0} while tt corresponding to time, XX to HBrO2,Y\mathrm{HBrO}_{2}, Y to Br−,Z\mathrm{Br}^{-}, Z to Ce4+\mathrm{Ce}^{4+}, and VV to BrMA. The rate constants and 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 kf[s−1]k_{f}\left[s^{-1}\right]. As certain values of this parameter can give rise to chaos, the following analysis was performed in order to identify different dynamics.

Table 2. Parameters of the investigated system (2).

4. Main results

The system of differential equations 2 was solved numerically in Matlab [18] using the ode45o d e 45 solver. The simulations were done depending on free parameter kfk_{f} ranging from 3×10−43 \times 10^{-4} to 5×10−45 \times 10^{-4} with 10−710^{-7} step. Each simulation was performed for the final time τ=100\tau=100 with a time step 10−410^{-4}. To avoid the systems distortions, only the last 20%20 \% of simulations were used for further computations. In all cases, initial conditions were set

x0=z0=v0=1x_{0}=z_{0}=v_{0}=1

As a main results, phase diagrams, amplitude frequency spectrums (FFT), and Poincaré sections were done for relevant choices of the parameter kfk_{f}. To illustrate changes in dynamical behavior, bifurcation diagrams underlined by the approximate entropy and the 0−10-1 test for chaos with suitable magnifications to the parameter kfk_{f} were plotted.

Consequently, as a goal of this paper, bifurcation diagrams, the approximate entropy, and the 0-1 test for chaos were computed for nested set of parameters kfk_{f}. The 0−10-1 test for chaos splits the values of the parameter for which regular (periodic or quasi-periodic) and irregular (chaotic) movements appear, while the output of approximate entropy detects increasing complexity of the investigated system 2 .
4.1. Phase diagrams, the Fourier spectra, and bifurcation diagrams. Periodic as well as chaotic dynamics was identified in the studied model (2). For example, in Fig. 1 and Fig. 2, regular movement is shown by a trivial loop (Fig. 1) for kf=3×10−4k_{f}=3 \times 10^{-4}, and a non-trivial loop (Fig. 2) for kf=3.2×10−4k_{f}=3.2 \times 10^{-4}. Figure 3 gives an example of a chaotic trajectory, that is for kf=3.5×10−4k_{f}=3.5 \times 10^{-4}.

The Fourier spectra were computed by Fast Fourier transform for kf=3×10−4,kf=3.2×10−4k_{f}=3 \times 10^{-4}, k_{f}=3.2 \times 10^{-4}, and kf=3.5×10−4k_{f}=3.5 \times 10^{-4}, shown in Figs. 1, 2, and 3 respectively. Regular behavior is observable for the first 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 number of harmonic frequencies, hence the frequency of the periodic trajectory is computable. Periodic motions of trajectory is also visible in Poincaré sections.

In the case of chaotic movement, Fig. 3, the Fourier spectra is formed by a number of harmonic components having the basic, super-harmonic, sub-harmonic, and combination frequencies on which further 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 is underlined by the Poincaré section.

Next, the bifurcation diagram (constructed as a projection of a local maxima) of the model (2) was plotted for each variable x,zx, z and vv with respect to the free parameter kf∈(3×10−4,5×10−4)k_{f} \in\left(3 \times 10^{-4}, 5 \times 10^{-4}\right)

Figure 1. Phase diagram with Poincaré section by plane v=0.8434v=0.8434 and Fourier transform of variable xx for kf=3×10−4k_{f}=3 \times 10^{-4}.

Figure 2. Phase diagram with Poincaré section by plane v=0.8445v=0.8445 and Fourier transform of variable xx for kf=3.2×10−4k_{f}=3.2 \times 10^{-4}.

Figure 3. Phase diagram with Poincaré section by plane v=0.847v=0.847 and Fourier transform of variable xx for kf=3.5×10−4k_{f}=3.5 \times 10^{-4}.
in Figure 4. In this bifurcation diagram, so-called “period doubling” and “windows” effects are also visible. Periodic movement can be identified in the range of the parameter, e.g., kf∈(3×10−4,3.24×k_{f} \in\left(3 \times 10^{-4}, 3.24 \times\right. 10−4)\left.10^{-4}\right) and kf∈(3.95×10−4,5×10−4)k_{f} \in\left(3.95 \times 10^{-4}, 5 \times 10^{-4}\right). The interval in between these values is interrupted by chaotic cases around kf=3.25k_{f}=3.25, and some kf∈(3.34×10−4,3.65×10−4)k_{f} \in\left(3.34 \times 10^{-4}, 3.65 \times 10^{-4}\right) and kf∈(3.85×10−4,3.9×10−4)k_{f} \in\left(3.85 \times 10^{-4}, 3.9 \times 10^{-4}\right).

Figure 4. Bifurcation diagram of the variable xx (left), zz (center), vv (right) kf∈k_{f} \in [3×10−4,5×10−4]\left[3 \times 10^{-4}, 5 \times 10^{-4}\right].
4.2. Approximate entropy. The approximate entropy is a technique used to quantify the amount of regularity and unpredictability of fluctuations in time-series. The main advantages are that it can be computed on short time series and it allows to compare the differences in complexity of the same system with different parameters settings, see, e.g., [21] to be detected. More complex notion of entropy type can be found in, e.g., [14]. To compute the approximate entropy, two

parameters must be set: embedding dimension mm and neighborhood threshold rr. Let s(t)∈Rs(t) \in \mathbb{R} for t={1,2,…,N}t=\{1,2, \ldots, N\} be a time series with NN observations. Then embedded vector S(t)S(t) at time tt, is defined as S(t)=[s(t),s(t+1),s(t+2),…,s(t+(m−1))]S(t)=[s(t), s(t+1), s(t+2), \ldots, s(t+(m-1))], where tt is the observed time and mm is the embedding dimension. The maximum distance of embedded vectors is computed as follows:

D(i,j)=d(S(i),S(j))=max⁡k=0,1,…,m−1∣s(i+k)−s(j+k)∣D(i, j)=d(S(i), S(j))=\max _{k=0,1, \ldots, m-1}|s(i+k)-s(j+k)|

for i,j={1,2,…,N−(m−1)}i, j=\{1,2, \ldots, N-(m-1)\}.
Compute the thresholded version of the distance with threshold given by rr :

dr(i,j)={1,D(i,j)<r0, otherwise d_{r}(i, j)= \begin{cases}1, & D(i, j)<r \\ 0, & \text { otherwise }\end{cases}

for i,j∈{1,2,…,N−(m−1)}i, j \in\{1,2, \ldots, N-(m-1)\}.
Compute Cim(r)C_{i}^{m}(r) as a ratio between points in the neighborhood of i and the number of embedded vectors.

Cim(r)=∑j=1N−(m−1)dr(i,j)N−(m−1)C_{i}^{m}(r)=\frac{\sum_{j=1}^{N-(m-1)} d_{r}(i, j)}{N-(m-1)}

Then compute the average of logarithm of all the Cim(r)C_{i}^{m}(r)

Φm(r)=1N−(m−1)∑i=1N−(m−1)ln⁡Cim(r)\Phi^{m}(r)=\frac{1}{N-(m-1)} \sum_{i=1}^{N-(m-1)} \ln C_{i}^{m}(r)

Finally, approximate entropy for the finite time series with NN data points is computed as

ApEn⁡(m,r,N)=Φm(r)−Φm+1(r)\operatorname{ApEn}(m, r, N)=\Phi^{m}(r)-\Phi^{m+1}(r)

For robust estimation, it was suggested by Pincus [21] the time series is containing at least 10310^{3} observations.

The approximate entropy was calculated using the TSEntropies package [28] for R [27]. The computations were made for the input vector ss given in a normalized form of all state variables:

s(t)=x2(t)+z2(t)+v2(t)s(t)=\sqrt{x^{2}(t)+z^{2}(t)+v^{2}(t)}

kf∈(3×10−4,5×10−4)k_{f} \in\left(3 \times 10^{-4}, 5 \times 10^{-4}\right) and r=0.1r=0.1. The results of approximate entropy for all values of the parameter kfk_{f} are in Figs. 7-9.
4.3. 0-1 test for chaos. The 0−10-1 test for chaos, invented by Gottwald and Melbourne [7], is one of the methods for distinguishing between regular and chaotic dynamics of a deterministic system. In contrast to the other approaches, the nature of the system is irrelevant, thus the test can be applied directly 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 to chaotic dynamics. With its easy implementation, evaluation, and wide range of application, using this tool for detecting chaos is becoming more popular in different fields.

The 0−10-1 test for chaos can be computed by the following algorithm [7].
Given the observation ϕ(j)\phi(j) for j=1,2,…,Nj=1,2, \ldots, N and a suitable choice of c∈(0,2π)c \in(0,2 \pi), the following translation variables are computed:

pc(n)=∑j=1nϕ(j)cos⁡(jc)qc(n)=∑j=1nϕ(j)sin⁡(jc)\begin{aligned} p_{c}(n) & =\sum_{j=1}^{n} \phi(j) \cos (j c) \\ q_{c}(n) & =\sum_{j=1}^{n} \phi(j) \sin (j c) \end{aligned}

Figure 5. The plot of pp versus qq for c=1.569853c=1.569853 for kf=3×10−4k_{f}=3 \times 10^{-4} showing regular dynamics (left) and chaotic dynamics for kf=3.5×10−4k_{f}=3.5 \times 10^{-4} (right).
for n=1,2,…,Nn=1,2, \ldots, N. The dynamics of the translation components pcp_{c} and qcq_{c} is shown on the pcp_{c} versus qcq_{c} plot. A bounded trajectory is in Fig. 5 (left) corresponding to regular movement, for kf=3×10−4k_{f}=3 \times 10^{-4}. An unbounded trajectory is in Fig. 5 (right) related to the chaotic case, for kf=3.5×10−4k_{f}=3.5 \times 10^{-4}.

The idea for the 0−10-1 test, first described in [7], is that the boundedness or unboundedness of the trajectory {(pj,qj)j∈[1,N]}\left\{\left(p_{j}, q_{j}\right)_{j \in[1, N]}\right\} can be studied through the asymptotic growth rate of its time-averaged mean square displacement (MSD), which is defined as

M(n)=lim⁡N→∞1N∑j=1Nd(j,n)2M(n)=\lim _{N \rightarrow \infty} \frac{1}{N} \sum_{j=1}^{N} d(j, n)^{2}

where

d(j,n)=(pj+n−pj)2+(qj+n−qj)2d(j, n)=\sqrt{\left(p_{j+n}-p_{j}\right)^{2}+\left(q_{j+n}-q_{j}\right)^{2}}

is the time lapse of the duration n(n≪N)n(n \ll N) starting from the position at time jj. As it is shown in [7,8][7,8], it is important to use values of nn small enough compared to NN, noted ncut,(n≤ncut)n_{c u t},\left(n \leq n_{c u t}\right). A subset of time lags ncut∈[1,N/10]n_{c u t} \in[1, N / 10] is advised for the computation of each KcK_{c}.

For bounded trajectories and regular dynamics, M(n)M(n) is a bounded function in time, whereas unbounded trajectories, meaning chaotic dynamics, are described by M(n)M(n) growing linearly with time. Thus the asymptotic growth rate of the MSD must be calculated, which correlates with the unboundedness of the trajectory.

As proposed in [7], the modified MSD is calculated as

D(n)=M(n)−E(ϕ)21−cos⁡(nc)1−cos⁡cD(n)=M(n)-E(\phi)^{2} \frac{1-\cos (n c)}{1-\cos c}

The output of the 0−10-1 test for chaos is computed by the correlation method as

Kc=corr⁡(ξ,Δ)∈[−1,1]K_{c}=\operatorname{corr}(\xi, \Delta) \in[-1,1]

for the vectors ξ=(1,2,…,ncut)\xi=\left(1,2, \ldots, n_{c u t}\right) and Δ=(Dc(1),Dc(2),…,Dc(ncut))\Delta=\left(D_{c}(1), D_{c}(2), \ldots, D_{c}\left(n_{c u t}\right)\right).
The final result of the test is

K=median⁡(Kc)K=\operatorname{median}\left(K_{c}\right)

Figure 6. The plot of KcK_{c} based on cc for kf=3×10−4k_{f}=3 \times 10^{-4} showing regular dynamics (left) and for kf=3.5×10−4k_{f}=3.5 \times 10^{-4} showing chaotic dynamics (right).

Figure 7. The result of approximate entropy (in blue) and the result of the 0-1 test for chaos (in red) for kf∈(3×10−4,5×10−4)k_{f} \in\left(3 \times 10^{-4}, 5 \times 10^{-4}\right). The magnification in the black rectangle is shown in Fig 8. The bifurcation diagram for variable xx is shown in the background.

The position of the studied system (2) at any moment of time is determined by displacements xx, zz, and vv, which are used for defining vector ϕ\phi :

ϕ(j)=x(j)2+z(j)2+v(j)2\phi(j)=\sqrt{x(j)^{2}+z(j)^{2}+v(j)^{2}}

For these simulations, a free software environment R [27] was used including the Chaos01 package developed by T. Martinović [17]. Comparison of values KcK_{c} for periodic and chaotic case is shown in Fig. 6, for kf=3×10−4k_{f}=3 \times 10^{-4} and kf=3.5×10−4k_{f}=3.5 \times 10^{-4}, respectively.

The results of the 0-1 test for chaos for all values of the parameter kfk_{f} are in Figs. 7-9.

Figure 8. The result of approximate entropy (in blue) and the result of the 0−10-1 test for chaos (in red) for kf∈(3.25×10−4,3.35×10−4)k_{f} \in\left(3.25 \times 10^{-4}, 3.35 \times 10^{-4}\right). The magnification in the black rectangle is shown in Fig 9. The bifurcation diagram for variable xx is shown in the background.

Figure 9. The result of approximate entropy (in blue) and the result of the 0−10-1 test for chaos (in red) for kf∈(3.322×10−4,3.324×10−4)k_{f} \in\left(3.322 \times 10^{-4}, 3.324 \times 10^{-4}\right). The bifurcation diagram for variable xx is shown in the background.

5. CONCLUSIONS

In this paper, the GF model (2) associated with the BZ chemical reaction was solved using adaptive six-stage, fifth-order, Runge-Kutta method implemented as ode 45 solver in Matlab. To eliminate the stiffness problem the model (2) was also simulated by ode 23s23 s solver in Matlab [18], outputs were identical.

The simulations were used to plot 3D phase portraits, bifurcation diagrams, the approximate entropy, and the 0−10-1 test for chaos. The mining process of dynamical properties were performed in the free software R [27] using packages TSEntopies [28] and Chaos01 [17] depending on flow rate parameter kfk_{f}.

It is evident from main results in Figs. 7-9 that both tests clearly detect regular and irregular patterns for given kfk_{f}.

Our results show that the method of approximate entropy returns qualification constant which describes complexity in the system invariantly with respect to the origin. On the other hand, the 0−10-1 test as qualification tool returns zero for regular (periodic or quasi-periodic) movement and one for irregular (chaos) characteristics.

Moreover, if the output of the 0−10-1 test is not close to zero or one, then the examined test case has not 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−10-1 test for chaos. In general, the increasing values of the 0−10-1 test for chaos are coupled to increasing approximate entropy and vice versa.

We notice isolated low values of the 0−10-1 test for chaos accompanied by comparatively low values of approximate entropy well within of chaotic region characterized by high 0−10-1 test chaos values and approximate entropy. To investigate and zoom in, we have constructed three-stage system of nested subintervals of flow rates kfk_{f}, see Figs. 7-9, for which in every level the 0−10-1 test for chaos and approximate entropy was computed. At every level we observe the same pattern.

That naturally yields suggestion of a fractal structure in the set of kfk_{f} :
Open Problem 1. Is there a totally disconnected (Cantor) set of flow rates kfk_{f} in [3×10−4,5×10−4]\left[3 \times 10^{-4}, 5 \times 10^{-4}\right] 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 kfk_{f} in [3×10−4,5×10−4]\left[3 \times 10^{-4}, 5 \times 10^{-4}\right] such that for each such parameter the GFG F model (2) is showing regular movement?

ACKNOWLEDGEMENT

This work was supported by The Ministry of Education, Youth and Sports from the National Programme of Sustainability (NPU II) project “IT4Innovations excellence in science - LQ1602”; by The Ministry of Education, Youth and Sports from the Large Infrastructures for Research, Experimental Development and Innovations project “IT4Innovations National Supercomputing Center - LM2015070”; by SGC grant No. SP2019/125 “Qualification and quantification tools application to dynamical systems”, VŠB - Technical University of Ostrava, Czech Republic, Grant of SGS No. SP2019/84, VŠB - 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. Patermaster, F. Rossi, J. Math. Chem. 56, 2876 (2018)
[3] I.R. Epstein, J.A. Pojman, An Introduction to Nonlinear Chemical Dynamics: Oscillations, Waves, Patterns, and Chaos (Oxford University Press, 1998)
[4] R.J. Field, E. Körös, 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örgyi, R.J. Field, J. Phys. Chem. 95, 6594 (1991)
[10] L. Györgyi, R.J. Field, Nature 355, 808 (1992)
[11] L. Györgyi, 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č, J. Math. Chem. 57, 1670 (2019)
[17] T. Martinovič, Chaos01: 0-1 Test for Chaos. (2016), https://CRAN.R-project.org/package=Chaos01\. Accessed 1 October 2018
[18] MATLAB 2015b (The MathWorks, Inc., 2015)

[19] O.V. Noskov, A.D. Karavsev, 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 Computing, 2018)
[28] J. Tomčala, 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ŠB - Technical University of Ostrava, 17. listopadu 15/2172, 70833 Ostrava, Czech Republic, Department of Applied Mathematics, VŠB - Technical University of Ostrava, 17. listopadu 15/2172, 70833 Ostrava, Czech Republic
E-mail address: judita.nagyova.st@vsb.cz
IT4Innovations, VŠB - Technical University of Ostrava, 17. listopadu 15/2172, 70833 Ostrava, Czech Republic
E-mail address: branislav.janzik@vsb.cz
IT4Innovations, VŠB - Technical University of Ostrava, 17. listopadu 15/2172, 70833 Ostrava, Czech Republic, Department of Applied Mathematics, VŠB - Technical University of Ostrava, 17. listopadu 15/2172, 70833 Ostrava, Czech Republic
E-mail address: marek.lampart@vsb.cz