Stability analysis for Selkov-Schnakenberg reaction-diffusion system

This study provides semi-analytical solutions to the Selkov-Schnakenberg reaction-diffusion system. The Galerkin method is applied to approximate the system of partial differential equations by a system of ordinary differential equations. The steady states of the system and the limit cycle solutions are delineated using bifurcation diagram analysis. The influence of the diffusion coefficients as they change, on the system stability is examined. Near the Hopf bifurcation point, the asymptotic analysis is developed for the oscillatory solution. The semi-analytical model solutions agree accurately with the numerical results.


Introduction
The emergence of oscillatory and multiple steady state solutions and chaotic behaviors is an interesting phenomenon observed in many chemical, biological, and physical systems. A number of experimental studies of chemical systems have been done to examine oscillating solutions involving Briggs-Rausher [1], Belousov-Zhabotinsky [2], and Bray-Liebhafsky [3] reactions, which demonstrate periodic variations in concentrations as visually striking color changes [4]. The Selkov-Schnakenberg system is an example of an oscillating systems associated with cellular processes in biochemical system [5]. The oscillatory phenomena in these systems have been investigated using continuous-flow well-stirred tank reactors (CSTRs), governed by ordinary differential equations (ODEs) for concentrations. CSTRs have exhibited excellent results in theoretical and experimental studies of chemical oscillations. Another reactor type important for the study of oscillating reactions is the reaction-diffusion cell. It is described using partial differential equations (PDEs) that can be used to characterize the emergence of many phenomena.
In 1979, Schnakenberg introduced a simple chemical reaction model for glycolysis that showed limit cycle behavior. The reaction scheme, known as Selkov-Schnakenberg model, occurs in the three steps: In the above reactions, A and B are two chemical sources, while U and V are autocatalyst and reactant, respectively. The most common examples of autocatalytic reactions are the chloride-iodide-malonic acid reaction and the reaction of phosphofructokinase glycolysis that includes adenosine triphosphate (ATP), adenosine diphosphate (ADP), and adenosine monophosphate (AMP) [6,7].
A great literature has been dedicated to the study of (1), including the Selkov model and the Schnakenberg model; see, for instance, [8][9][10][11][12][13][14] and the references therein. A generalized Selkov-Schnakenberg reaction-diffusion system is analyzed in [15]. The authors provide values for the stability of the unique constant steady state. Also, the effect of the diffusion coefficients on the stability of the system has been investigated in this study. In [16], the Selkov-Schnakenberg model is considered, and the steady state problem is studied. They obtained results show a change in the stability of the system, leading to the formation of different spatial patterns. In [17], the Selkov-Schnakenberg system was studied numerically, and numerical results for snaking of patterns were obtained.
One of the most important mathematical methods that has proven effective in developing semi-analytical models for reaction-diffusion systems is the Galerkin method [18,19]. It is based on using the approximation of a system of PDEs with a system of ODEs with the same behavior. This method was used in many research studies, for example, in 2002, Marchant [20] employed this method to analyze Gray-Scott model. In this study, the regions of stability and instability in the system were determined using combustion theory. This method yielded higher accuracy compared to the semi-analytical model and the PDE system (its numerical results). Also, in [21,22] this method was applied to the reversible Selkov model. The researchers found the Hopf regions and bifurcation diagram for the model, as well they studied the effect of feedback control in the boundary conditions of the system with and without precursor and final products. Moreover, [23,24] also considered a semi-analytical model for the Schnakenberg reaction-diffusion system. The two studies focused on the stability and singularity behavior of the model, showing that the model has a hysteresis bifurcation pattern. To further elaborate on the application of the Galerkin method in many reaction-diffusion problems, we recommend the reader to refer to these references [25][26][27][28][29][30][31][32].
The purpose of this study is to investigate the singularity behavior and stability of the Selkov-Schnakenberg system which, through previous studies, have not yet been investigated and so all the results of this work are novel and genuine. This paper is structured as follows. Section 2 shall present the mathematical model of the Selkov-Schnakenberg system and describes the methodology for application of the Galerken method to PDEs and boundary conditions. Section 3 uses steady state concentration profiles and bifurcation diagrams to illustrate the complexity of the steady state multiplicity of the system. Section 4 reports the results of the application of the singularity theory to determine the hysteresis bifurcation points and defines the region of the two bifurcation diagram patterns. In 5, stability analysis of the semi-analytical model is performed and Hopf bifurcation region is determined. In Section 6, an asymptotic analysis of the periodic solution near the Hopf bifurcation point in the semi-analytical model is performed.
2 The semi-analytical model

The mathematical model
The Selkov-Schnakenberg system (1) is governed by the following coupled nonlinear PDEs in a one-dimensional ( D 1 − ) reaction-diffusion cell: In this system, u u x t , = ( ) and v v x t , = ( ) are the dimensionless concentrations of the autocatalyst u and reactant v, respectively. The diffusion coefficients for u and v are represented by the parameters D u and D v , respectively. While δ 2 is a positive constant, λ and δ 1 are non-negative constants. The reactor has permeable boundaries, x 1 = , linked to a reservoir containing both u and v in zero concentrations, thus the system is considered open. At x 0 = , the center of the cell, the concentration flux across the boundary condition is zero. When the two parameters λ and δ 1 are zeros, the model (2) converts to the Selkov model, suggested by Selkov in 1969 [33,34], while if λ 0 = and δ 0 1 > , the model becomes the Schankenberg model [35]. The numerical solution of (2) is obtained utilizing the Crank-Nicolson method for solving the PDEs numerically with the accuracy of O t x Δ , Δ 2 ( ).

The Galerkin method
The semi-analytical model of the Selkov-Schnakenberg system is obtained through the use of the Galerkin method in a 1-D reaction-diffusion cell. This converts the PDEs (2) governing the system into a system of ODEs. A function expansion c o s c o s , is used, where θ πx = . This expansion fulfills the spatial boundary conditions (3). At the same time, it satisfies the PDEs (2) in an average sense, but not exactly. Moreover, the expansion used in (4) has the property that at the cell center the concentrations are u u u 1 2 = + and v v v 1 2 = + . We define the parameters shown in the expansion equations by evaluating averaged versions of PDEs system and then multiplying by the basis function. As a result, the following system of ODEs is produced: The system of ODEs is found by truncating the series of basis functions (4) after two terms. This truncation is sufficient to show high accuracy in the results without enlarging the expression. To ensure sufficient accuracy, it is compared to the one-term solution obtained by assuming u 0 2 = and v 0 2 = in (4).

Analysis of steady state bifurcations
This section deals with the steady state bifurcations of the system (5). To obtain the steady state solutions, the time derivatives  It may be found that the solution for the two-term model is more accurate than the one for the one-term model, using the numerical solution as a reference, with relative error rates of less than 0.43 and 0.65% for u and v, respectively. , and δ 0.7 2 = . The solutions for both models (oneand two-term) together with the numerical solution of PDEs system are drawn and displayed. For these parameter values, the curves demonstrate a unique pattern, and over the whole range of concentrations, only a single steady state solution may be found for the system. It may be found from the figures that u increases and v decreases as λ increases. On the other hand, the results in Figure 3 It may be found that once again the results of the two-term model and the numerical solution to the PDEs agree to a great extent.

Application of singularity theory
Singularity theory is based on a theoretical study that describes cases of steady state behavior in the systems of ODEs [36]. Several studies have examined the feasibility of applying this theory to chemical reactions. One of the most famous of these studies is [37] in which the conditions for designating regions of isola and hysteresis bifurcation curves in CSTR were presented.
This section aims to apply the theory of singularity to the semi-analytical model (2) to detect the steady state bifurcation diagrams that can emerge in this system. The two primary bifurcation types are those with unique and breaking wave patterns. It is possible to write the steady state equations for the two-term semianalytical model in the general form (4) as follows: Hence, our selected bifurcation parameter is λ. The unique and breaking wave solutions illustrated in the previous section correspond to the emergence of the hysteresis bifurcation curve. The defining conditions for finding the hysteresis singularity points, described in details by [20], for the two-term semianalytical model are: In calculations, we regard u v v , , 2 1 2 , and λ as implicit functions of u 1 . The hysteresis singularity conditions for the one-term model are obtained as follows: Figure 4 shows the hysteresis bifurcation curve between the autocatalyst diffusion coefficient D u and the reactant diffusion coefficient D v , which divides the plane into two different areas. The figure shows the hysteresis points for both one-and two-term models, which are obtained by solving the conditions (7) and (8) the figure that increasing the autocatalyst diffusion coefficient D u causes a shift from the breaking wave to unique pattern for fixed values the reactant diffusion coefficient D v values. The example in Figure 2, parameter values D 0.5 u = and D 0.04 v = , lies above the relevant curve shown in Figure 4 and hence is a unique solution, while the example shown in Figure 3, parameter values D 0.1 u = and D 0.05 v = , is located below the curve and hence is a breaking wave solution. Figure 5 demonstrates the two-term hysteresis curve with δ 0.005, 0.03 1 = , and 0.08. It shows that increasing the rate constant δ 1 leads to the hysteresis curve moving upward in the D D u v − parameter space. In brief, the altering in the rate of the constant δ 1 leads to the appearance or disappearance of multiplicity of steady state.
, and all the other parameters are given in Figure 2.  and δ 0.7 2 = .

The phenomena of limit cycles and Hopf bifurcations
The limit cycles and Hopf bifurcations are phenomena considered very important for studying physical, biological, and chemical systems. The Hopf bifurcation theory was explained in a variety of sources on dynamical systems and bifurcations theory [38,39]. The Hopf bifurcation results from the emergence of a periodic solution due to local switching of the properties of the steady state solutions from stable to unstable when a simple complex conjugate pair of eigenvalues crosses over the complex plane's imaginary axis. This section explores the local stability of the semi-analytical model to investigate the effect of the diffusion coefficients on the system stability (2). The degenerate Hopf points are computed to predict a semianalytical map in which Hopf bifurcations and limit cycles occur. A comparison of this forecast with numerical results is then made. First, the stability of the ODEs system (5) representing the two-term semi-analytical model is studied. The Jacobian matrix is the following: When one pair of the eigenvalues is completely imaginary, Hopf bifurcations are observed: I a a a a a a 0.
Therefore, solving the following condition yields the degenerate Hopf points: Results in Figure 6 show the degenerate Hopf bifurcation curve for the two-term model of (11)  , and δ 0.08 1 = . The parameter is fixed as δ 0.7 2 = .
region of the graph, stable solutions exist, and there is no limit cycle. In the lower portion of the plot, the system is destabilized; unstable solutions can occur and the limit cycle arises. It may be found in Figure 6 that the increase in the diffusion coefficient of the autocatalyst concentration D u corresponds to a decrease in the diffusion coefficient of the reactant concentration D v , which leads to expansion of the stable region and contraction of the Hopf bifurcation region. For example, at D 0.04 v = , the solution for the two-term semi-analytical model is D 0.292 u = , whereas the numerical solution is D 0.293 u = ; therefore, the error amount is less than 0.4%. Hence, it is seen that the change in diffusion coefficients affects the stability of the system. = − + . In the two previous figures, the results of the two-term model yield high accuracy as referred to the numerical solution. Figure 8(a) shows that the limit cycle occurs in a phase plane of the autocatalyst concentration u against the concentration of the reactant v. Figure 8 curve in Figure 6, where the compatibility of the emergence of the limit cycles with the theoretical prediction in the ODEs system is shown. The period of oscillation in the numerical solution is estimated as 9.3, while it is equal to 9 and 9.4 for one, and two-term semi-analytical models, respectively. Moreover, the twoterm semi-analytical model amplitudes of oscillation are 1.    Figure 8. These results demonstrate that the decreasing λ causes the expansion of the limit cycle, while increasing λ causes the lessening of the amplitude of the limit cycle, leading to a steady state scenario. Figure 11 shows the phase plane u v − representing the two-term limit cycle curves for different values of parameter δ 1 = 0.004, 0.06 and 0.1. Other parameters are fixed as λ 0.06 , and δ 0.7 2 = . The figure shows that the increase in the rate constant δ 1 causes a reduction in the amplitude of the limit cycle. In  summary, different values of both parameters λ or δ can lead either to stability or to perturbation of the system.

Oscillatory solutions near the Hopf bifurcation point
This section elaborates on the periodic solutions of the semi-analytical ODEs model (5) for the Selkov-Schnakenberg system and applies asymptotic analysis to them. The Lindstedt-Poincaré method [40] perturbation theory is used to compute asymptotic solutions to both the one-and two-term semi-analytical models.
The determination of periodic solutions of power series with respect to a small oscillation amplitude parameter is considered. The method is applied to remove secular terms and obtain corrections to the bifurcation parameter and frequency [41].

Lindstedt-Poincaré perturbation method for the one-term model
One-term semi-analytical model comprises the following system of ODEs: The Lindstedt-Poincaré method is utilized to find the π 2 -periodic solution by expanding the solution of the system (12) in the following form: where τ ωt = is the dimensionless time, ω is the frequency of the oscillator, and ε is known as the parameter of periodic oscillations amplitude which can be found by the following normalization condition:  , and 0.1. Other parameters are as in Figure 8.
Then, the frequency ω and the bifurcation parameter λ are expanded in a power series of ε 2 as λ λ ε λ ω ω ε ω , , where λ 0 and ω 0 are the parameter values at the Hopf bifurcation point. The perturbation corrections λ 2 and ω 2 can be computed by utilizing the solvability conditions in the third order of the perturbation analysis. Inserting the expansion equations (13) and (15) into (12) yields the following perturbation equations for the first three orders of ε  ( ) For the first two orders, the solutions of (16) and (17) where c.c. is the complex conjugate. By substituting the solutions (19) into (16) and (17), and then splitting into real and imaginary parts, we can find the amplitudes A A iA is not a root of homogeneous equations. Consequently, the right value of C along with λ 2 and ω 2 will cancel out secular terms in O ε 3 ( ).
Hence, the limit cycle as well as its period may be expressed as where the extrema of the oscillatory solutions are given by u ε 0.693 2 1 = ± and v ε 1.686 3.822 1 = ± . The limit cycle (22) demonstrates the classical square root behavior near the Hopf Point λ 0.119 = .
In addition, the oscillation frequency increases linearly with the increase in λ. A Hopf bifurcation point takes place at λ 0.119 = .
6.2 Lindstedt-Poincaré perturbation method for the two-term model The ODEs system for the two-term semi-analytical model is found in (5). The periodic solutions near the Hopf bifurcation point are obtained through the same procedure as for the one-term model, but with some additional complexities. The Lindstedt-Poincaré solution of (5) has the following form: , .
The perturbation equations for the first three orders of ε are found by substituting (23) and (15) Again, substituting (24) and (25)   Stability analysis for Selkov-Schnakenberg reaction-diffusion system  59 The leading-order periodic solution is  (27) where the extrema of the periodic solutions are u ε 0.677 2.054 = ± and v ε 1.269 4.022 = ± . In this case, the Hopf bifurcation point occurs at λ 0.104 = . Figure 12 demonstrates the bifurcation diagram for the concentration of autocatalyst u and the reactant v for the one-term solution. Figure 13 illustrates the two-term solution in the same way. The semi-analytical and perturbation results in the two figures are shown for the same parameter values: D 0.09 u = , D 0.02 v = , δ 0.03 1 = , and δ 0.7 2 = . It may be inferred from the bifurcation diagram that the system has a stable state. This state is attained after the Hopf bifurcation point: λ 0.119 = for the one-term and λ 0.104 = for the twoterm. This response can be expected based on the bifurcation diagram illustrated in Figures 12 and 13, indicating that the system has oscillatory solutions for a sufficiently small rate of λ; here, the system initially shows a limit cycle. Shortly after, the system crosses the instability threshold resulting in the suppression of the limit cycle. Both figures show a reliable comparison of the results of the semi-analytical and Lindstedt-Poincaré perturbation methods for λ 0.103 > . The validity range, however, is limited due to the rapid growth of oscillatory amplitude and the limited accuracy range for the expansion of the Taylor series, with respect to decreasing λ.

Conclusion
Within the scope of this paper, a detailed study of the Selkov-Schnakenberg system in the reaction-diffusion cell is presented. The PDEs system of the model was approximated with the system of ODEs using the Galerkin method. The steady state and bifurcation diagrams for the system have been constructed based on the singularity theory and discussed. The Hopf bifurcation region was also studied and the limit cycles in the system have been shown. An asymptotic analysis of the periodic solution near the Hopf bifurcation point was carried out and its results were compared to the semi-analytical solution. It has been shown that the diffusion coefficients for both autocatalyst and reactant, as well as the λ and δ 1 parameters, all affect the stability of the system. Unparalleled compatibility of the semi-analytical and numerical solutions was shown throughout the examples in this paper.
This method has shown its effectiveness and accuracy and may be applied to other chemical models, studying the stability of which is associated with greater complexity such as the three-component FitzHugh-Nagumo model and the three-component reversible Gray-Scott model. The effect of feedback and time delay in boundary conditions on the stability of the considered model should be further investigated.