Stability analysis for Selkov-Schnakenberg reaction-diffusion system

1 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

3 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 d d t u i = d d t v i = 0 , i = 1 , 2 in versions of (5) were set to zero and the four simultaneous equations ( f i = 0 , i = 1 , … , 4 ) were solved numerically. For the one-term model, the two simultaneous equations were obtained by equating u 2 = v 2 = 0 .

Figure 1(a) and (b) depict the steady state concentration profiles for the autocatalyst and the reactant, u and v , against x . It shows the solutions for both models (one- and two-term) with the numerical solution of PDEs system. The parameters used here are λ = 0.01 , D u = 0.4 , D v = 0.07 , δ 1 = 0.03 , and δ 2 = 0.7 . For both the autocatalyst u and reactant v , concentrations of the solutions have a rounded peak inside the reactor. Note that the concentrations peak at the center of the cell, at x = 0 , where ( u , v ) = ( 0.0477 , 4.833 ) for the one-term model and ( u , v ) = ( 0.0463 , 4.643 ) for the two-term model, while ( u , v ) = ( 0.0465 , 4.673 ) for the numerical solution. 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.

Figure 1

Steady state concentration profiles against x for (a) autocatalyst u and (b) reactant v , respectively. All parameter values were fixed at λ = 0.01 , D u = 0.4 , D v = 0.07 , δ 1 = 0.03 , and δ 2 = 0.7 .

The bifurcation diagrams demonstrate the steady state concentrations of both u and v in the cell center at x = 0 versus the bifurcation parameter λ . In Figure 2(a) and (b), steady state bifurcation diagrams of u and v concentrations vs λ are shown. The values used for the other parameters here are D u = 0.5 , D v = 0.04 , δ 1 = 0.03 , and δ 2 = 0.7 . The solutions for both models (one- and 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(a) and (b) show the pattern of the breaking wave with the same parameters in Figure 2(a) and (b) except for D u and D v which are 0.1 and 0.05, respectively. The steady state location folds back to itself when λ is in the domain (0.003, 0.008). Between the two turning points, the system has a multiplicity of steady state solutions where the autocatalyst concentration suddenly increases. At the same time, the reactant concentration drops at the bifurcation point λ = 0.007 . 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.

Figure 2

The unique pattern of steady state bifurcation diagrams for (a) autocatalyst u and (b) reactant v concentrations, respectively, against the bifurcation parameter λ , at x = 0 . All parameter values were fixed at D u = 0.5 , D v = 0.04 , δ 1 = 0.03 , and δ 2 = 0.7 .

Figure 3

The breaking wave pattern of steady state bifurcation diagrams for (a) autocatalyst u and (b) reactant v concentrations, respectively, against the bifurcation parameter λ , at x = 0 . In these graphs, D u = 0.1 , D v = 0.05 , and all the other parameters are given in Figure 2.

4 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 semi-analytical 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 semi-analytical model are:

In calculations, we regard u 2 , v 1 , v 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), with δ 1 = 0.03 and δ 2 = 0.7 . The two regions of the graph represent two different bifurcation diagrams. The lower portion of the plot has the breaking wave pattern with a region of multiple steady state solutions. At the same time, in the upper region of the plot exists a unique steady state solution. It may be found from 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 u = 0.5 and D v = 0.04 , 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 u = 0.1 and D v = 0.05 , is located below the curve and hence is a breaking wave solution.

Figure 4

The D u − D v parameter plane shows the degenerate hysteresis curve which separates the domains of different steady state patterns. The parameters are fixed as δ 1 = 0.03 and δ 2 = 0.7 .

Figure 5 demonstrates the two-term hysteresis curve with δ 1 = 0.005 , 0.03 , and 0.08. It shows that increasing the rate constant δ 1 leads to the hysteresis curve moving upward in the D u − D 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.

Figure 5

The D u − D v parameter plane shows the two-term hysteresis curve with δ 1 = 0.005 , 0.03 , and δ 1 = 0.08 . The parameter is fixed as δ 2 = 0.7 .

5 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 semi-analytical 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:

The eigenvalues, Γ , of J are computed using the characteristic equation and growth of small perturbations in the steady state solution. The following quartic equation is the characteristic equation:

When one pair of the eigenvalues is completely imaginary, Hopf bifurcations are observed:

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 the numerical results of governing PDEs in the plane of the autocatalyst diffusion coefficient D u versus the reactant diffusion coefficient D v . The other pair of parameter values for the emergence of the Hopf bifurcation is δ 1 = 0.03 and δ 2 = 0.7 . The degenerate Hopf curve separates the plane into two regions. In the upper 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 v = 0.04 , the solution for the two-term semi-analytical model is D u = 0.292 , whereas the numerical solution is D u = 0.293 ; 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.

Figure 6

The locus of the degenerate Hopf curve in the D u − D v parameter plane is shown. Other parameters remain as in Figure 4.

Figure 7 shows another possible example. It displays the loci of Hopf bifurcation lines for both the two-term model and the numerical solution in the ( δ 1 , δ 2 ) plane for D u = 0.4 and D v = 0.03 . Again, the degenerate Hopf line divides the domain into two regions. In the domain above the line, there is a steady state solution and no Hopf points, while in the domain below the line, solutions are unstable and the limit cycle occurs. Note that, the Hopf points occur below the line so that δ 1 ∈ ( 0 , 0.16 ) and δ 2 ∈ ( 0.55 , 0.73 ) for physical solutions. In summary, it can be concluded that an increase in the concentration of the constant δ 1 makes the degenerate Hopf bifurcation line move down towards the smaller values of the constant δ 2 . This means the stability region of the system expands. Moreover, the Hopf bifurcation line may be expressed with the following linear equation: δ 2 = − 1.025 δ 1 + 0.721 . In the two previous figures, the results of the two-term model yield high accuracy as referred to the numerical solution.

Figure 7

The δ 1 − δ 2 parameter plane representing the locus of the degenerate Hopf curve at D u = 0.4 and D v = 0.03 .

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(b) and (c) illustrate the change in time of the autocatalyst concentration u and the reactant concentration v , respectively. It may be found from these figures that the system stabilizes to its oscillatory dynamics, allowing estimating the period of the solution and amplitudes of the concentrations. The parameter values utilized are λ = 0.06 , D u = 0.1 , D v = 0.01 , δ 1 = 0.03 , and δ 2 = 0.7 . Figure 8 displays the results of both one- and two-term models along with the numerical solution. The values of the parameters used in this example lie in the domain below the Hopf bifurcation 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 two-term semi-analytical model amplitudes of oscillation are 1.74 and 2.68 for the concentrations of autocatalyst and reactant, respectively, whereas the numerical amplitudes of concentrations are 1.70 and 2.65, respectively. At the same time, using the parameter values λ = 0.02 , D u = 0.5 , D v = 0.06 , δ 1 = 0.03 , and δ 2 = 0.7 , lying in the domain under the Hopf bifurcation curve in Figure 6, results in a steady state solution. These time series of both autocatalyst and reactant concentrations, respectively, are shown in Figure 9(a) and (b). When time exceeds 40, the system reaches the steady state which is given by ( u , v ) ≃ ( 0.073 , 5.181 ) for one-term model, ( u , v ) ≃ ( 0.072 , 4.961 ) for the two-term model, and ( u , v ) ≃ ( 0.072 , 4.995 ) for the numerical solution. The results presented in Figures 8 and 9 present further evidence for the quality and accuracy of the two-term method, with the numerical solution as a reference point.

Figure 8

The limit cycle curve within the u against v phase plane (a) and the time evolutions for u (b) and v (c), respectively, at x = 0 . The parameter space values used are λ = 0.06 , D u = 0.1 , D v = 0.01 , δ 1 = 0.03 , and δ 2 = 0.7 .

Figure 9

The time evolutions for u (a) and v (b) against t at x = 0 . The parameter values used are λ = 0.02 , D u = 0.5 , D v = 0.06 , δ 1 = 0.03 .

Figure 10 shows phase portraits of the concentrations at the center of the reactor for various choices of bifurcation parameter λ , illustrating the limit cycle curves. The figure presents the results of the two-term model at λ = 0.06 (the blue closed curve), λ = 0.08 (the black closed curve), λ = 0.1 (the red closed curve), and λ = 0.115 (the green closed curve), while other parameter values are the same as in 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 , D u = 0.1 , D v = 0.01 , and δ 2 = 0.7 . 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.

Figure 10

Phase plane u − v illustrating the limit cycle curves for different values of bifurcation parameter λ = 0.06 , 0.08 , 0.1 , and 0.115. Other parameters are as in Figure 8.

Figure 11

Phase plane u − v illustrating the limit cycle curves for different values of δ 1 = 0.004 , 0.06 , and 0.1. Other parameters are as in Figure 8.

6 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].

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:

(23) u 1 ( τ ) = u 1 s + ε u 11 ( τ ) + ε 2 u 12 ( τ ) + ε 3 u 13 ( τ ) , v 1 ( τ ) = v 1 s + ε v 11 ( τ ) + ε 2 v 12 ( τ ) + ε 3 v 13 ( τ ) , u 2 ( τ ) = u 2 s + ε u 21 ( τ ) + ε 2 u 22 ( τ ) + ε 3 u 23 ( τ ) , v 2 ( τ ) = v 2 s + ε v 21 ( τ ) + ε 2 v 22 ( τ ) + ε 3 v 23 ( τ ) .

The perturbation equations for the first three orders of ε are found by substituting (23) and (15) into (5). By doing so, four equations for each order of ε are obtained (due to their excessively long expressions, they are not shown here). The form of the solution for the first order is

(24) u 11 ( τ ) = e i τ + c.c. , v 11 ( s ) = A 1 e i τ + c.c. , u 21 ( τ ) = A 2 e i τ + c.c. , v 21 ( τ ) = A 3 e i τ + c.c.

To find the complex amplitudes A 1 , A 2 , and A 3 , (24) is substituted into the O ( ε ) equations. The solution for the order two O ( ε 2 ) is

(25) u 12 ( τ ) = B 1 e 2 i τ + c.c. , v 12 ( τ ) = B 2 e 2 i τ + c.c. , u 22 ( τ ) = B 3 e 2 i τ + c.c. , v 22 ( τ ) = B 4 e 2 i τ + c.c.

Again, substituting (24) and (25) into the O ( ε 2 ) equations results in four equations for complex amplitudes B i . In the third order equations, secular terms, e ± i τ , appear in the right-hand side. To cancel them out, we set ( u 13 , v 13 , u 23 , v 23 ) = C 1 e i τ ( 1 , 1 , 1 , 0 ) T and choose C 1 , λ 2 , and ω 2 appropriately.

As a result, the solution for the special case is obtained; it is expounded below. The Hopf bifurcation arises at

(26) ( λ 0 , ω 0 , u 1 s , v 1 s , u 2 s , v 2 s ) = ( 0.104 , 0.726 , 0.683 , 1.926 , − 0.006 , − 0.657 ) , A 1 = − 1.109 + 1.607 i , A 2 = 0.025 + 0.072 i , A 3 = − 0.248 − 0.123 i , B 1 = − 0.047 + 1.472 , B 2 = − 1.475 + 1.234 i , B 3 = − 0.009 + 0.322 i , B 4 = − 0.505 − 0.495 i , α 2 = − 2.361 , ω 2 = − 2.144 .

The leading-order periodic solution is

(27) u ( τ ) ≃ u 1 s + u 2 s + ε ( 2 cos ( ω t ) + ℜ { A 2 e i τ + A ¯ 2 e − i τ } ) , v ( τ ) ≃ v 1 s + v 2 s + ε ( ℜ { A 1 e i τ + A ¯ 1 e − i τ + A 3 e i τ + A ¯ 3 e − i τ } ) , ε ≃ 0.044 − 0.424 λ , ω ≃ 0.909 λ + 0.632 ,

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 u = 0.09 , D v = 0.02 , δ 1 = 0.03 , and δ 2 = 0.7 . 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 two-term. 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 λ .

Figure 12

Bifurcation diagrams for (a) autocatalyst u and (b) reactant v concentrations, respectively, against bifurcation parameter λ . The one-term semi-analytical and perturbation results for the parameter values D u = 0.09 , D v = 0.02 , δ 1 = 0.03 , and δ 2 = 0.7 are shown.

Figure 13

Bifurcation diagrams for (a) autocatalyst u and (b) reactant v concentrations, respectively, against bifurcation parameter λ . The two-term semi-analytical and perturbation results for the parameter values D u = 0.09 , D v = 0.02 , δ 1 = 0.03 , and δ 2 = 0.7 are shown.

7 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.