Dynamics of the Stochastic Belousov-Zhabotinskii Chemical Reaction Model

: In this paper, we discuss the dynamic behavior of the stochastic Belousov-Zhabotinskii chemical reaction model. First, the existence and uniqueness of the stochastic model’s positive solution is proved. Then we show the stochastic Belousov-Zhabotinskii system has ergodicity and a stationary distribution. Finally, we present some simulations to illustrate our theoretical results. We note that the unique equilibrium of the original ordinary differential equation model is globally asymptotically stable under appropriate conditions of the parameter value f , while the stochastic model is ergodic regardless of the value of f .


Introduction
The theoretical analysis of repeated oscillation processes in open systems began in Lotkal [1,2]. Bak [3] and Higgins [4] considered both closed and open systems and generalized the theory of such reactions and Spangler and Snell [5] considered a model system. Interest in oscillating reactions was generated by the large number of such processes observed in biological systems and these ideas were developed by Prigogine and coworkers [6,7]. Oscillating reactions generally involve two autocatalytic reactions that turn each other off and on and are somewhat analogous to the famous flip-flop circuit in electronics. Field et al. [8] showed that it is possible to devise a mechanism of oscillation involving only a single autocatalytic process.
The iodine catalyzed decomposition of hydrogen peroxide [9] was the first reaction believed to oscillate homogenously. Belousov observed the second example of oscillating chemical reactions in homogeneous solutions, i.e., oxidation of tetravalent-trivalent cerium ion coupled catalytic citric acid by potassium bromate. The Belousov-Zhabotinskii reaction was found in 1950s by the former Soviet Union biophysics Belousov and Zhabotinskii [10]. So we called the reaction as Belousov-Zhabotinskii reaction, referred to as B-Z reaction. In 1969, Prigogine proposed the theory of dissipative structure, which clearly explains the reason for the occurrence of oscillation reaction. This makes the B-Z reaction back to the focus of the study. The theory of dissipative structure argues that when the system is far from equilibrium state, that is, it is in the state of nonlinear and non-equilibrium, the disordered homogeneous state is not necessarily stable.
There are many explanations for the mechanism of B-Z reaction [11,12]. In one version, the total reaction consists of two processes: A and B. The total reaction of processe A is: The total reaction of processe B is: Process A is the ion reaction of double electron transfer, process B involves free radicals and single electron transfer. When the concentration of Bromide ion exceeds the critical concentration process A arises mainly, and when the concentration of Bromide ion below the critical concentration process B arises mainly. Bromide ion is the control species here. In the process, B bromate ion not only oxidate the metal ions, but also generate bromine. Such oscillation can sustain thousands of times in the closed system, and the reaction process does not need to add the reactants. So this kind of reaction provides great convenience for the study of the chemical wave. The reaction solutions appear in two colors. If the solution is shallow, it will be similar to the phenomenon of wave interference, two kind colors of the wave spread alternately. The wave in the same color disappeared when the contact happened, which is different from electromagnetic wave.
In 1990, L.Gyorgy, T.Turanyi, R.J. Field [8] elaborated a skeletal reaction mechanism. The reactions constituting the mechanism of the oscillatory B-Z reaction may be divided into an inorganic and an organic subset. It contains eighty elementary reactions, and it is so complex that peoples put forward a simplified mechanism of six steps, which is expressed as: The reaction substrates included bromate, cerium ammonium sulfate (or cerium sulfate), malonic acid and dilute sulfuric acid, and bromate is a substrate which can not be changed. Metal ions is generally the Ce or Mn, but can also be the complex ions formed by Fe, Ru, Co, Cu, Cr, Ag, Ni, Os. Malonic acid can be instead of other reducer. The formation of other halogen-containing oxysalt and anion, such as chloride ions, can interfere with the reaction. B-Z reaction usually join the adjacent ferroin as the indicator, it is the complexes of phenanthroline and ferrous ion, which is red in the reduced state and blue in the oxidation state. Cerium ion is yellow in oxidation state and trivalent cerium ion is colorless in the reduced state. So the effect of synthesis is: the oxidation state green, reducing state red. Now we will introduce the establishment of the mathematical model for the B-Z chemical reaction and consider the concentration changes of some main reactants. First, we will explain why we have not considered the acidity or the water concentration in this paper? It is because their variations in the media are negligible due to the use of some tampon. , Q is the concentration of waste generation. We describe the above reactions in the following five steps: (i) Describe the transform from Y to X; (ii) Inactivating X and Y at the same time; (iii) and (iv) autocatalytic product X; (v) Bimolecular decomposition of X; (vi) The comprehensive reaction change Z to Y. We express these steps as: where f is a stoichiometric factor of an overall reaction step, and k i (i = 1, · · · , 6) is the reaction rate constant which contains the H + concentration effect.
From the basic principle of modelling, we get the mathematical model of the reaction which is often referred to as B-Z reaction model: Note as in [13] with the dimensionless transformation system (1) is transformed to: Tyson et al. [14,15] studied the dynamics of model (3) and system (3) has two equilibria O(0, 0, 0) and P(x * , y * , z * ), and O is unstable. Here, The globally asymptotic stability of the equilibrium state can be seen and in [13,16]. The Theorem 1.10 in [13] shows that the positive equilibrium state of model (2) is globally asymptotic stable for f ∈ (0, f 1 ) ( f 2 , ∞). In [17], Tang proved that for system (3), B is a positive invariant set, and all positive initial solution will enter into B ultimately, here Note model (1) is deterministic and does not incorporate the effect of environment noise, which is always present. In fact, chemical reaction models are affected by environmental white noise, because the temperature and pressure varies during chemical reaction processes and is closely related to the process. Two stochastic chemical reaction models were discussed in [18][19][20] where the chemical reaction model was considered under stochastic perturbation with the Lyapunov method.
In this paper, we introduce stochastic perturbations into the system (1) and we consider the stochastic system: where W 1 (t), W 2 (t), W 3 (t) are independent Brownian motions and δ 2 1 > 0, δ 2 2 > 0, δ 2 3 > 0 represent the intensities of white noise. Consider the dimensionless transformation of this stochastic model and the relation of the coefficient as in (2) except that Then system (4) becomes: where The aim of this paper is to study the dynamical behavior of system (5). In Section 2, we show the solution of system is positive and global. In Section 3, we prove that system (5) has ergodicity and a stationary distribution and the last section presents some simulations to illustrate our theoretical results.

Existence and Uniqueness of the Positive Solution
It is necessary that the solutions of the stochastic model is positive in a chemical reaction model, so we prove the existence and uniqueness of the positive solution for model (5).
Unless otherwise specified, we let (Ω, F , {F t } t≥0 , P) be a complete probability space, with a filtration {F t } t≥0 satisfying the usual conditions (i.e., it is right continuous and F 0 contains all P-null sets) throughout this paper. Denote Generally, we consider the following d-dimensional stochastic differential equation with initial value and B(t) denotes the standard d-dimensional Brownian motions defined on the above probability space. The differential operator L associated with the upper equation is defined by If the differential operator L acts on a function V ∈ C 2,1 (S h × R + ; R + ), then we have
Proof. In order to illustrate the positive local solution (x(τ), y(τ), z(τ)), τ ∈ [0, τ e ) of stochastic B-Z system (5) is global, we just need to prove that τ e = ∞ a.s. Let m 0 ≥ 0 be an enough large number such that both x 0 , y 0 , z 0 lie within the interval [ 1 m 0 , m 0 ]. Set m ≥ m 0 , and define the stopping time we set inf φ = ∞ (as usual φ denotes the empty set) throughout this paper. Clearly, τ m is increasing s. If we prove τ ∞ = ∞ a.s., then τ e = ∞ obviously, i.e., (x(0), y(0), z(0)) ∈ R 3 + a.s. for all t ≥ 0. Let me put it another way, tend to complete the proof we just need to prove that τ e = ∞ a.s. If not, there are T > 0 and ε ∈ (0, 1), such that Therefore, there exists an integer m 1 ≥ m 0 bring that P{τ m ≤ T} ε f or all m ≥ m 1 .
are two positive constants. The non-negativity of the C 2 −function V can be seen from u − 1 − log u ≥ 0, for all u > 0. Applying Itô's formula, we can derive that Choosing the value of k 1 , k 2 as in (7), then

Ergodicity
Before we start to prove the ergodicity of stochastic B-Z reaction model, a result which can be found in [21] (we also refer the reader to [22]) will be present.
Let X(t) be a homogeneous Markov Process, which have been described by the stochastic equation in E l (E l denotes the euclidean l-space) and its diffusion matrix is The differential operator L, which associated with Equation (9) is defined by

Lemma 1.
(See [21]) Suppose that there is a bounded domain U ⊂ E l with regular boundary Γ, which having the following characters: (i) The smallest eigenvalue of the diffusion matrix Λ(x) is bounded away from zero in the domain U and some neighborhood thereof.
(ii) If x ∈ E l \U, the mean time τ at which a path issuing from x reaches the set U is finite, and sup x∈K E x τ < ∞ for every compact subset K ⊂ E l . Then the Markov process X(t) has a stationary distribution µ(·). If we let f (·) be a function integrable with regard to the measure µ. Then for all x ∈ E l .

Remark 1.
The proof of Lemma 1 is given in [21]. In Theorem 4.  [23], Chapter 3, p. 103 and Rayleigh's principle in [24], Chapter 6, p. 349). In order to verify (ii), it is enough to show that there are some neighborhood U and a non-negative C 2 -function, therefore for any E l \U, LV is negative (for details see [25], p. 1163).

Lemma 2.
Let X(t) be a regular homogeneous temporally Markov process in E l . In case X(t) is recurrent relative to a few bounded domain U, then it is recurrent relative to any nonempty domain in E l .

Remark 2.
Due to the existence of a positive solution of stochastic B-Z model (5) was obtained in Theorem 1, it is sufficient to take R 3 + as the whole space.
Proof. In order to prove this theorem, all we need to do is to verify that (i) and (ii) hold according to Lemma 1. First, the stochastic system (5) can be written as: Here the diffusion matrix is which implies condition (i) is satisfied. Now we need to verify condition (ii). First, we construct a nonnegative C 2 −function V and a closed set U ∈ ∑ (which lies in R 3 + entirely) so as to sup (x,y)∈R 3 + \U LV(x, y) < 0, which guarantees that condition (ii) is satisfied. Take into account a C 2 −function h(x, y, z): h(x, y, z) = x + 2s 2 y + 3s f ω z − l 1 log x − l 2 log y − l 3 log z, (x, y, z) ∈ R 3 + .
Here l 1 , l 2 , l 3 are positive numbers with 0 < l 1 < 1 and Obviously, if σ 2 1 < 2s, we only need to choose l 2 , l 3 small enough such that condition (10) holds. It is not hard to verify that h(x, y, z) has a unique minimum point (l 1 , 2s 2 . Next we define a nonnegative C 2 −function taking the following form: ).
By direct calculation, we obtain that 2 ), and from (11) we have C < 0. Then we can get Define a closed set and ε is a sufficiently small number such that − s(1 − l 1 ) 1 − s f 1 Let For any (x, y, z) ∈ R 3 + \U we talk over the following six cases: {− 1 2 qsx 2 + (s + 3s f + l 1 sq + l 2 s )x}. Then we can get LV < −1 on D ∞ 1 from (12) and (13).
Case 6. On D 3 , using the same method as in Case 5, we have the following inequality: which gives LV < −1 in this domain, since (12) and (18) hold. Thus condition (ii) is satisfied. Then we finish the proof of Theorem 1.

Simulation
Now, we test our theory conclusions by simulations in this section. In this section we suppose that the unit of time is a minute and reactant concentration are measured in units of mol/L·min. We always use a discretization method and choose t = 0.002. Choose the parameters in the system (5) as respectively. Obviously, only the value of the parameter f is different in the above two conditions. We set the initial value (x(0), y(0), z(0)) = (0.8, 0.8, 0.5) and simulate the scatter distribution and the stationary distribution figures of the stochastic differential equation model (5) with Matlab under the condition (19) and (20) . The simulation results are shown in Figures 1-4 respectively. The corresponding deterministic model (3), the value of the parameters in (19) do not satisfy the condition in Theorem 1.10 [13], because the value f = 1 satisfies neither (1) nor (2). The value of f = 0.1 in (20) satisfies the conditions of Theorem 1.10 [13]. In these two cases, the stochastic model (5) always has ergodicity and a stationary distribution. In this paper, we have added multiplicative noise terms into the Equation (1) as a representation of environmental white noise, furthermore we proved that the existence of the positive solution and the stochastic model is ergodic. In fact, ergodicity is a kind of weak stability, that is, stability according to distribution. Regardless of what state the system show, the chemical reaction is always continuing in the weak stable state.    Figure 3. The solution of the stochastic system and its histogram when parameters are chosen as in (19). The blue lines represent the solution of the corresponding undisturbed system (3) and the red lines represent the solution of the stochastic system (5). The right picture is the histograms of the stochastic system (5).  (20). The blue lines represent the solution of the corresponding undisturbed system (3) and the red lines represent the solution of the stochastic system (5). The right picture are the histograms of the stochastic system (5).

Conclusions
The conclusion we get is very interesting that original restrictions on parameter f in the deterministic is gone, all of the condition is the white noise satisfies inequality σ 2 1 < 2s. So we think white noise is conducive to the stability of the B-Z reaction stochastic system in practice. In future work, we will try our best to make a breakthrough on this model.