Turing Instability of Brusselator in the Reaction-Diffusion Network

Turing instability constitutes a universal paradigm for the spontaneous generation of spatially organized patterns, especially in a chemical reaction. In this paper, we investigated the pattern dynamics of Brusselator from the view of complex networks and considered the interaction between diffusion and reaction in the random network. After a detailed theoretical analysis, we obtained the approximate instability region about the diffusion coefficient and the connection probability of the random network. In the meantime, we also obtained the critical condition of Turing instability in the network-organized system and found that how the network connection probability and diffusion coefficient affect the reaction-diffusion system of the Brusselator model. In the end, the reason for arising of Turing instability in the Brusselator with the random network was explained. Numerical simulation verified the theoretical results.


Introduction
Pattern formation, a kind of nonuniform macroscopic structure with some regularity in space or time, is ubiquitous. It was first proposed by Turing [1] for systems containing morphogens, although initially, they may be very uniform, a pattern or structure may later emerge due to the instability of the uniform equilibrium, which is triggered by random disturbances. A theoretical analysis of Turing instability in a semidiscrete Brusselator model has been investigated [2]. In [3], the authors proposed and discussed the Brusselator model in a random framework. e mean-field equation, which proved that an organized Turing pattern could be produced in a specific parameter region, was derived. To reveal the effect of external noise on the system, a detailed random analysis of the Brusselator scheme was conducted. e stochastic analysis revealed that such systems' structural stability would be disturbed even if the bifurcation parameters are subject to small external disturbances. ere would be different space and time structures in a certain range of noise intensity and correlation time [4]. It is well known that the Brusselator model is a typical model for studying patterns, so it has attracted the interest of many scholars in different fields. Firstly, the superdiffusion term on pattern formation and pattern selection in the Brusselator model has been well studied. ey found that Turing instability can occur under superdiffusion even though the diffusing initiator is faster than the inhibitor [5]. And the model with superdiffusion has been well investigated in [6]. For the Brusselator model, the nonlinear diffusion term and the linear diffusion on the pattern patterns are compared. e process of pattern formation in one-dimensional and two-dimensional space domain was also studied in [7]. In addition, the Brusselator model also attracts interest from scholars studying chemistry and stochastic oscillations in [8,9]. en, from a basic point of equilibrium, the authors studied the regular Hopf bifurcation and the singular Hopf bifurcation of the Brusselator model under the reaction periodic force and obtained a method suitable for the study of the nonlinear vibration of periodic forces in general [10]. Besides, the existence of the pattern formation of the Brusselator model under homogeneous Neumann boundary conditions and the prediction of pattern formation caused by Turing instability under certain assumptions have been well studied [11]. e study of random networks has become more and more popular among scholars in recent years. It is widely used in various fields such as biology, chemistry, and engineering, as discussed in [12][13][14]. In the study of chemical reactions, random networks can be used to describe the self-diffusion of molecules. In [15], the relationship between the eigenvalues of the Laplace matrix and the degree was pointed out. e localization properties of the Laplace eigenvectors in various random networks were explained. In particular, it provided a theoretical basis and method for further study of the dynamic characteristics of correlated stochastic networks. In [16], McCullen and Wagenknecht emphasized the importance of network structure, revealed the basic connection between small-scale activity patterns on the network and local pattern formation of the whole discipline, and investigated the reaction-diffusion system on the complex network topology. Regarding the analysis of self-organized systems on the network, Nakao and Mikhailov provided a new perspective on analyzing the pattern formation of activator-inhibitor systems on the network in [17]. e system instability caused by the diffusion term on an undirected random network with a certain probability was studied. e pattern theory of the directed network can be derived from [18]. en, the Turing instability of the reaction-diffusion model defined on the complex network was studied in [19], and three types of models on the complex network were shown in the article. Numerical results showed that the uniform steady-state stability region depends on the network system's structure in the diffusion coefficient space. In [20], they inferred a general theoretical proof that the Turing system's three key features are directly determined by the topology. Recently, Mimar et al. [21] proved that the degree of Laplace can reflect the system's topological characteristics under different networks, which is related to the local characteristics of the diffusion coefficient.
Nowadays, Zheng and Shen [22] investigated the pattern formation in FitzHugh-Nagumo model with a random network, obtained the approximated Turing instability region about the diffusion coefficient and connection probability, and gave a feasible method for studying reactiondiffusion system with connection probability. Although it is known that connection probability plays an essential role in random networks, the influence of random networks on Brusselator cannot be ignored. Still, there are a few literature studies on the effects of network connection probability on Brusselator. In this paper, we explored the effect of random networks on the pattern formation of the reaction-diffusion system from node connection probability.
Next, we will combine the above methods to investigate the pattern dynamics behaviour of the Brusselator model in a random network with connection probability. For the network Laplacian matrix, the connection probability between network node pairs also plays an extremely significant role. Besides, the connection probability affects the diffusion term of the system by changing Laplacian eigenvalues and then affects the stability of the reaction-diffusion system. In Section 2, based on the positive equilibrium point, the Brusselator model's stability is analyzed. e critical condition of Turing instability for the reaction-diffusion system concerning the diffusion coefficient was obtained. In Section 3, the theoretical analysis of Turing instability for the reactiondiffusion system with diffusion term was introduced. In Section 4, the theoretical results obtained in the paper are summarized, and the results are verified by numerical simulation.

Linear Stability Analysis of Brusselator.
e reaction between molecules is described by the Brusselator model, which is a mathematical model proposed by the Brussels school to simulate self-organized phenomena. e Brusselator model studied in this paper is given by the following reaction formula [4]: where k i (i � 1, 2, 3, 4) is a positive parameter representing the reaction rate constant. According to the law of mass action, the differential equation of X and Y concentration can be written as follows: where c M is the constant variable representing the concentrations of M (M is the substance A, B.), and c X and c Y are the variables representing the concentrations of X and Y.
To dimensionless, the differential equation (2) becomes the following: where x, y, t, a, and b are scale variables, and their expressions are as follows: It is well known that for chemical reaction systems, the positive equilibrium point has scientific significance. erefore, we have a great deal of interest in the nonnegative equilibrium point and focus on the system's stability near the positive constant fixed points. Obviously, system (3) has the unique constant solution M(x, y), where x � a and y � (b/a). According to the coordinate translation transformation, the equilibrium point M is translated to the origin O(0, 0), that is, introducing X � x + x, Y � y + y in the system. en, the linearized system (3) can be expressed as follows: To analyze the stability of the equilibrium point of system (5), we write the corresponding Jacobian matrix of the linearized system as follows: e characteristic equation of the system at point O is given by According to Weida's theorem and local equilibrium point stability theory, the equilibrium O is stable if and only if P < 0 and Q > 0 holds. From Q � a 2 and the parameter domain, it is obvious that Q > 0, and from p � b − 1 − a 2 , we can obtain parameters satisfying the condition which is consistent with the condition that system (3) has a unique positive equilibrium point.

Lemma 1. If the condition (9) holds, the unique positive equilibrium point M(x, y) of system (3) is asymptotically stable.
Proof. From the stability analysis of the above linearized system, it is clear that when condition (9) holds, the real part of the eigenvalues is negative, that is, the equilibrium point O is stable. erefore, the equilibrium point M of nonlinear system (3) is the stable focus. When condition (9) holds, the unique positive constant solution M of the initial reaction system is asymptotically stable (Figure 1). e proof is completed. □

Brusselator Model with Random Network.
e Brusselator model with the diffusion term can be described as the following equation set: Many scholars have well studied the system with the nonnetwork diffusion term. In this paper, the random network can describe the self-diffusion of molecules to investigate the stability change of the homogeneous state of system (3), where the connecting probability of a pair of network nodes is p. e following steps generate the random network and the adjacency matrix element: Step 1: suppose the network consists of n nodes Step 2: the value of the element A ij of the adjacency matrix is generated as follows: there is an edge between i and j when the random number <p, which is Figure 2 shows a random network structure with a connection probability of p � 0.05. en, the balance equation (3) corresponding to each network node can be rewritten as follows: e general solution of equation (3) can be expressed as follows:

Complexity 3
Firstly, we consider the system's unique positive equilibrium state M(x, y) after joining the network diffusion, en the Jacobian matrix of each node in the network diffusion system is expressed as follows: en the characteristic function of system (11) can be written as follows: so the necessary and sufficient condition for Turing instability is that there exists at least In addition, in order to obtain the Turing instability region of the system, we analyzed the stability of the reaction-diffusion system and obtained its critical value . It is obvious that this equation has two roots k 2 1c , k 2 2c , when Reλ(k 2 ) � 0. When the continuous system is unstable, Reλ(k 2 ) � 0 is established in Figure 4, that is, the value range is In addition, the relationship between the network node degree k i and the local eigenvectors − Λ i of the Laplacian matrix and the eigenvalues can be obtained in Figure 5. e detailed proofs of Lemmas 2 and 3 are given in reference [22].

Application of Mean-Field Approximation.
en, the mean-field approximation theory is applied to analyze further the reaction-diffusion system, which can explain the Turing instability mechanism induced by the network diffusion term. In the mean field, the reaction-diffusion system can be expressed as follows:   4 Complexity where C x � N j�1 A ij x j , C y � N j�1 A ij y j , and k i represents the degree of the network.
Mean-field approximation is a method to study complex multibody problems, turning multibody problems into monomer problems. For the Brusselator model, the interaction between a single molecule and other molecules is replaced by an external field effect on this molecule. erefore, when the role of other molecules is fixed, and only the equilibrium state M(x, y) of the system is considered, the corresponding a single node system can be expressed as follows: From the linear analysis of system (16), the characteristic equation is given as follows: From the above analysis of the stability of Brusselator, system (16) without the diffusion effect is stable when 0 < b < 1 + a 2 holds. Assuming that λ 1 and λ 2 are the eigenvalues of equation (17), the stability of system (16) depends on the sign of

Results
In this section, we present an explanation of theoretical results based on chemical mechanisms. Firstly, the network adjacency matrix A is a symmetric matrix which is randomly generated based on a random network with probability p. e relationship between the Laplacian matrix and the Complexity 5 adjacency matrix of the network is L ij � A ij − k i δ ij , where k i is the degree of the node i and δ ij is a Dirac function, if i � j, δ ij � 1, otherwise δ ij � 0. We choose the parameters a � 1 and b � 1.5 and derive the relationship between the diffusion coefficients d 1 and d 2 when the Hopf bifurcation occurs from Figure 3, which shows the critical condition for d 1 and d 2 when Hopf instability occurs. When diffusion coefficients d 1 and d 2 meet the critical value, Hopf instability occurs after control parameters are set. And we can get that the system may be in an instability state when d 2 > (10 + 4 � 6 √ )d 1 . When the ratio of d 2 to d 1 exceeds the critical value, the substance continues to interact to destroy the system's existing equilibrium state. However, when we consider the connection probability, not all values of the instability region can induce Turing instability in that the connection probability of network nodes will also affect the stability of the system. Next, we can give the numerical simulation to verify that Turing instability occurs which is not only related to the diffusion coefficients but also to the probability p. In addition, since the eigenvalues of the Laplacian matrix are discrete, it is a feasible method to study the effect of the distribution of eigenvalues on the system when we analyze continuous systems. From Lemma 2 and Figure 5, we can gain that the network Laplacian matrix eigenvalue − Λ i is proportional to the node degree so that it can be approximated by the degree of the node. e bifurcation diagram of the system concerning d 2 , displayed in Figure 6, shows that the system remains stable when p � 0.0004 and p � 0.69. And we perform some numerical simulations to certify the above analysis of the system's instability conditions, and the chemical mechanism of Turing instability in the system is given.
Firstly, we obtain the bifurcation diagram (Figures 7-11(a)) of the reaction-diffusion system with respect to d 2 when d 1 � 0.01 under different parameter p, which shows that the bifurcation point is consistent with the critical value for theoretical analysis.
In Figure 7, we can see that when p � 0.0005 > 0.0004, d 2 exceeds the critical value of instability; Turing instability will occur in the reaction-diffusion system. Figure 7(a) shows the bifurcation diagram of the system with respect to the diffusion coefficient, indicating that the system's stability will be destroyed with the increase of d 2 . Figure 7(b) shows that the system will appear unstable when d 2 � 1.5 and p � 0.0005. Figure 8 demonstrates the stability of the network diffusion system with respect to the diffusion coefficient at p � 0.1. Figure 8(a) shows that the bifurcation diagram of system (11) on the diffusion coefficient can be obtained when the connection probability is given. When d 2 � 0.1, Figure 8(b) demonstrates that the system is stable. But when p � 0.1 and d 2 � 0.3, the corresponding pattern formation (Figure 8(c)) is unstable. When p � 0.1 and d 2 � 0.3, the real part of the eigenvalues of the system characteristic equation (17) changes with k i and − Λ i is shown in Figure 8(d), which shows that the system is unstable and Turing instability occurs.
When p � 0.2, the bifurcation diagram of x i on d 2 is shown in Figure 9(a), which illustrates that the system's equilibrium point is stable when d 2 is lower than the critical bifurcation value. Still, Turing instability occurs when d 2 is larger than the critical bifurcation value, and the substance will continue to react. From the pattern formation, which is shown in Figures 9(b) and 9(c), it can be seen that the Brusselator system is stable when p � 0.2 and d 2 � 0.1 and the system is unstable when d 2 � 0.3, which is consistent with the bifurcation diagram on d 2 . From Figure 9(d), it can be concluded that when d 2 � 0.3, there exist some Re(Λ i ) falling into the unstable region, so the system is unstable.
When p � 0.33, the bifurcation diagram of x i on d 2 is shown in Figure 10(a), which illustrates that the system's equilibrium point is stable when d 2 is lower than the critical bifurcation value. Still, Turing instability occurs when d 2 is larger than the critical bifurcation value, and the substance will continue to react. From the pattern formation, which is shown in Figures 10(b) and 10(c), it can be seen that the Brusselator system is stable when p � 0.33 and d 2 � 0.1 and the system is unstable when d 2 � 0.3, which is consistent with the bifurcation diagram on d 2 . Figure 10(d) shows that when d 2 � 0.3, there are some Re(Λ i ) falling into the unstable region, so the system is unstable.
For p � 0.43, the bifurcation diagram of concentration x i with respect to d 2 can be drawn as Figure 11, and the bifurcation point is consistent with the Turing instability threshold. e bifurcation diagram shows that when d 2 is less than the critical value, the system is stable, otherwise unstable. As can be seen from Figures 11(b) and 11(c), the corresponding pattern information is stable when d 2 � 0.1, but unstable when d 2 � 0.3. In Figure 11(d), there exist some − Λ i belonging to D, so Turing instability occurs when When p � 0.53, we can give the bifurcation diagram of d 2 , as shown in Figure 12(a), and the pattern formation when d 2 � 0.1 is shown in Figure 12(b), which shows the system is stable. When the diffusion coefficient increases to a critical value, the uniform equilibrium point of the system begins to become unstable, and the current equilibrium state of the system is broken. So from Figure 12(c), when d 2 � 0.3, Turing instability occurs, and the system continues to interact. In Figure 12(d), there are some − Λ i falling into the instability region, so the system appears the Turing bifurcation phenomenon when d 2 � 0.3.
But when p � 0.8, the bifurcation of d 2 and the corresponding pattern formation of the system are exhibited in  Figure 13(d) demonstrates that the system's eigenvalues do not fall in the Turing instability region Λ ∩ D � Φ, so when d 2 � 0.3, the system is stable. As the connection probability increases, the system's bifurcation threshold will be increased or even disappeared. From Figure 4, we can conclude that Turing instability occurs when k i ∈ [8, 39] by the mean-field theory. From Figure 14(a), we can obtain that when k i � 7 < [8, 39] does not belong to the instability region, the system will remain in a stable state. When k i � 20 (or 30) ∈ [8, 39] falls into the unstable region, the stability of the equilibrium changes. Figures 14(b) and 14(c) mean that the molecular concentration does not reach a steady state, and the reaction continues. Figure 14(d) shows that the system is stable when k i � 50 > [8,39].
In Figure 15, we give the bifurcation diagram of the deterministic system (16) with respect to the node degree. e figure shows that the system has a transcritical bifurcation about the parameter k i near the equilibrium point,     where the bifurcation points are, respectively, k 1c � 7.033 and k 2c � 37.86, further confirming the system instability region about k i . For how the connection probability p affects the stability of the Brusselator system, the numerical simulation shows that the Turing instability occurs when the connection probability is very large or very small. e approximate instability region of the reaction-diffusion system for the Brusselator model with respect to d 2 and p is shown in Figure 16. erefore, from Lemma 3 and the above analysis when there is some Re(λ) in the Turing instability region Λ ∩ D ≠ Φ, Turing instability occurs. It means that the introduction of a random network disrupts the equilibrium state of the reaction system.
Finally, through the above analysis and numerical simulation, we obtained the approximate region of Turing instability of d 2 and p, which is given by Figure 16. In addition, we can also draw that the approximate instability region of the system on connection probability is (1/n 2 ) < p < (d 2c /d 1 )(4/5)(ln n/n), where (1/n 2 ) means that the network node pairs are almost sparse. Here, we only give the approximate instability region for the connection probability. However, the more precise instability region and the derivation process for the Brusselator reaction-diffusion system still need further study.

Conclusion
In conclusion, firstly, the Turing instability critical value of the Brusselator model with a random network is obtained in Figure 3. In addition, for the connection probability, through numerical simulation and the eigenvalue characteristics of the network matrix, we can conclude that the critical value of Turing instability of the system about p is either very large or very small, that is, the approximate range of p is (1/n 2 ) < p < (d 2c /d 1 )(4/5)(ln n/n). When the control parameters and d 1 are given, the bifurcation of the reactant concentration on d 2 under different connection probabilities verifies the approximate region of the system about p. It can be seen from the bifurcation graph that Turing instability will occur when p falls into the instability region from Figures 7-13(a). And we give the numerical simulation of the relation between the real part of the root of the characteristic equation and the eigenvalues of the network matrix. e   Figure 16: e instability region is related to the parameters d 2 and p, where the red region is the instability region. stability of the system near the equilibrium point is judged and verified by the real part of the eigenvalue exceeding the zero line in Figures 8-13(d) [18]. However, the upper bound of Turing instability on the connection probability is the best result of numerical simulation in this system, and its exact value needs to be explicitly derived. erefore, we need to study the precise value of the upper bound on the Turing instability of the system concerning p.
Finally, we apply the mean-field theory to transform the problem of multimolecule interactions into a single-molecule problem. And the instability region of the reactiondiffusion system with respect to the network node degree is obtained. When k i � 7 < [8, 39], the system remains stable at the equilibrium point. When k i � 20(30) ∈ [8, 39], the system has Turing instability. When k i � 50 > [8, 39], the system is stable at the positive equilibrium point. And the bifurcation of k i Figure 15 further verified this conclusion. In future research, we can further study the effect of random networks with specific characteristics on the reaction-diffusion system and work to determine the appropriate random networks for each system and the differences and connections between the effects of different random networks on the same reaction system.

Data Availability
e data used to support the findings of the study are available from the corresponding author upon request.

Conflicts of Interest
e authors declare that there are no conflicts of interest.