# Turing Instability of Brusselator in the Reaction-Diffusion Network

Authors: Yansu Ji; Jianwei Shen
Journal: Complexity
 (2020)

---

## Abstract

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.


---

## Body

## 1. 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. The 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. The stochastic analysis revealed that such systems’ structural stability would be disturbed even if the bifurcation parameters are subject to small external disturbances. There 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. They 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. The 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]. Then, 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].The 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–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. The 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]. The system instability caused by the diffusion term on an undirected random network with a certain probability was studied. The pattern theory of the directed network can be derived from [18]. Then, 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 reaction-diffusion 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 Section2, based on the positive equilibrium point, the Brusselator model’s stability is analyzed. The 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 reaction-diffusion 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.


## 2. The Analysis of the Brusselator Model

### 2.1. Linear Stability Analysis of Brusselator

The reaction between molecules is described by the Brusselator model, which is a mathematical model proposed by the Brussels school to simulate self-organized phenomena. The Brusselator model studied in this paper is given by the following reaction formula [4]:(1)A⟶k1X,B+X⟶k2Y+D,2X+Y⟶k33X,X⟶k4E,where ki (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:(2)dcXdτ=k1cA−k2cB+k4cX+k3cX2cY,dcYdτ=k2cBcX−k3cX2cY,where cM is the constant variable representing the concentrations of M (M is the substance A,B.), and cX and cY are the variables representing the concentrations of X and Y. To dimensionless, the differential equation (2) becomes the following:(3)dxdt=a−b+1x+x2y=fx,y,dydt=bx−x2y=gx,y,where x, y, t, a, and b are scale variables, and their expressions are as follows:(4)t=k4τ,a=k12k3k43cA,b=k2k4cB,x=k3k4cX,y=k3k4cY.It is well known that for chemical reaction systems, the positive equilibrium point has scientific significance. Therefore, 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 Mx¯,y¯, where x¯=a and y¯=b/a. According to the coordinate translation transformation, the equilibrium point M is translated to the origin O0,0, that is, introducing X=x+x¯,Y=y+y¯ in the system. Then, the linearized system (3) can be expressed as follows:(5)dXdt=b−1X+a2Y,dYdt=−bX−a2Y.To analyze the stability of the equilibrium point of system (5), we write the corresponding Jacobian matrix of the linearized system as follows:(6)J=b−1a2−b−a2.The characteristic equation of the system at pointO is given by(7)λi2−Pλi+Q=0,(8)P=b−1−a2,Q=a2.According to Weida’s theorem and local equilibrium point stability theory, the equilibriumO is stable if and only if P<0 and Q>0 holds. From Q=a2 and the parameter domain, it is obvious that Q>0, and from p=b−1−a2, we can obtain parameters satisfying the condition(9)0<b<1+a2,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 Mx¯,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. Therefore, 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). The proof is completed.Figure 1
When the parametersa=1 and b=1.5, the equilibrium of the Brusselator system (2) is asymptotically stable. (a) The time evolution diagram of the initial value around the equilibrium point. (b) The phase plane.
(a)(b)


### 2.2. Brusselator Model with Random Network

The Brusselator model with the diffusion term can be described as the following equation set:(10)∂xdt=fx,y+d1∇2x,∂ydt=gx,y+d2∇2y.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. The following steps generate the random network and the adjacency matrix element:Step 1: suppose the network consists ofn nodesStep 2: the value of the elementAij of the adjacency matrix is generated as follows: there is an edge between i and j when the random number <p, which is Aij=Aji=1, otherwise Aij=Aji=0. For example, Figure 2 shows a random network structure with a connection probability of p=0.05.Figure 2
Network structure is displayed when connection probabilityp=0.05.Then, the balance equation (3) corresponding to each network node can be rewritten as follows:(11)dxidt=fxi,yi+d1∑jLijxi,dyidt=gxi,yi+d2∑jLijyi,where Lij=Aij−kiδij is the Laplacian matrix, ki is the degree of note i, δij=0 if i≠j, and δij=1 otherwise.The general solution of equation (3) can be expressed as follows:(12)xi=∑m=1Ncmβmeλmtϕim,yi=∑m=1Ncmeλmtϕim,where ∑jLijϕjm=Λmϕim.Firstly, we consider the system’s unique positive equilibrium stateMx¯,y¯ after joining the network diffusion, where fx¯,y¯=gx¯,y¯=0. Suppose that the eigenvalues of the Laplacian matrix is Λii=1,2,…,n, which is arranged in descending order, i.e., 0=Λ1>Λ2>⋯>Λn. Then the Jacobian matrix of each node in the network diffusion system is expressed as follows:(13)Ji=b−1+d1Λia2−b−a2+d2Λi.Then the characteristic function of system (11) can be written as follows:(14)λ2−b−1−a2+d1Λi+d2Λiλ+d1d2Λi2+b−1d2−a2d1Λi+a2=0.From Section2.1, we can know that b−1−a2+d1Λi+d2Λi<0, so the necessary and sufficient condition for Turing instability is that there exists at least one Λi to satisfy wΛi=d1d2Λi2+b−1d2−a2d1Λi+a2<0 or Reλ>0. 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 satisfies 4d1d2a2−b−1d2−a2d12/4d1d2=0 (Figure 3). It is obvious that this equation has two roots k1c2,k2c2, when Reλk2=0. When the continuous system is unstable, Reλk2=0 is established in Figure 4, that is, the value range is k2∈D=λk1c2<λ<k2c2. In addition, the relationship between the network node degree ki and the local eigenvectors −Λi of the Laplacian matrix and the eigenvalues can be obtained in Figure 5.Figure 3
A region diagram of the occurrence of Turing instability for diffusion coefficient.Figure 4
Linear stability analysis. The relationship graph betweenReλ about Λik2 when a=1, b=1.5, p=0.1, and d1=0.01. The critical values are Λc=−2+6/10+46 and d2c=0.1+0.046.Figure 5
The role ofki. (a) The relationship graph between degree ki and eigenvalue −Λi. (b) About the dispersion relationship of ki is shown.
(a)(b)Lemma 2.
For the network Laplacian matrixL, generally kmax=maxki,kmin=minki and Λ is the eigenvalue of L. Therefore, Λ∈C=Λ−2kmax<Λ<M0−kmin.Lemma 3.
For a network-organized system, the system is always stable when all the eigenvalues of the network Laplacian matrix are not in the instability regionΛ∩D=Φ (empty set), the instability occurs when Λ∩D≠Φ, and D∩C≠Φ.The detailed proofs of Lemmas2 and 3 are given in reference [22].


### 2.3. Application of Mean-Field Approximation

Then, 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:(15)dxidt=fxi,yi+d1Cx−kixi,dyidt=gxi,yi+d2Cy−kiyi.where Cx=∑j=1NAijxj,Cy=∑j=1NAijyj, and ki 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. Therefore, when the role of other molecules is fixed, and only the equilibrium stateMx¯,y¯ of the system is considered, the corresponding a single node system can be expressed as follows:(16)dxidt=fxi,yi+d1x¯−kixi,dyidt=gxi,yi+d2y¯−kiyi.From the linear analysis of system (16), the characteristic equation is given as follows:(17)λ2−b−1−a2−d1ki−d2kiλ+d1d2ki2+a2d1−b−1d2ki+a2=0.From the above analysis of the stability of Brusselator, system (16) without the diffusion effect is stable when 0<b<1+a2 holds. Assuming that λ1andλ2 are the eigenvalues of equation (17), the stability of system (16) depends on the sign of wki, where wki=d1d2ki2+a2d1−b−1d2ki+a2. If wki>0, system (16) is stable, otherwise unstable.


## 2.1. Linear Stability Analysis of Brusselator

The reaction between molecules is described by the Brusselator model, which is a mathematical model proposed by the Brussels school to simulate self-organized phenomena. The Brusselator model studied in this paper is given by the following reaction formula [4]:(1)A⟶k1X,B+X⟶k2Y+D,2X+Y⟶k33X,X⟶k4E,where ki (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:(2)dcXdτ=k1cA−k2cB+k4cX+k3cX2cY,dcYdτ=k2cBcX−k3cX2cY,where cM is the constant variable representing the concentrations of M (M is the substance A,B.), and cX and cY are the variables representing the concentrations of X and Y. To dimensionless, the differential equation (2) becomes the following:(3)dxdt=a−b+1x+x2y=fx,y,dydt=bx−x2y=gx,y,where x, y, t, a, and b are scale variables, and their expressions are as follows:(4)t=k4τ,a=k12k3k43cA,b=k2k4cB,x=k3k4cX,y=k3k4cY.It is well known that for chemical reaction systems, the positive equilibrium point has scientific significance. Therefore, 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 Mx¯,y¯, where x¯=a and y¯=b/a. According to the coordinate translation transformation, the equilibrium point M is translated to the origin O0,0, that is, introducing X=x+x¯,Y=y+y¯ in the system. Then, the linearized system (3) can be expressed as follows:(5)dXdt=b−1X+a2Y,dYdt=−bX−a2Y.To analyze the stability of the equilibrium point of system (5), we write the corresponding Jacobian matrix of the linearized system as follows:(6)J=b−1a2−b−a2.The characteristic equation of the system at pointO is given by(7)λi2−Pλi+Q=0,(8)P=b−1−a2,Q=a2.According to Weida’s theorem and local equilibrium point stability theory, the equilibriumO is stable if and only if P<0 and Q>0 holds. From Q=a2 and the parameter domain, it is obvious that Q>0, and from p=b−1−a2, we can obtain parameters satisfying the condition(9)0<b<1+a2,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 Mx¯,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. Therefore, 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). The proof is completed.Figure 1
When the parametersa=1 and b=1.5, the equilibrium of the Brusselator system (2) is asymptotically stable. (a) The time evolution diagram of the initial value around the equilibrium point. (b) The phase plane.
(a)(b)


## 2.2. Brusselator Model with Random Network

The Brusselator model with the diffusion term can be described as the following equation set:(10)∂xdt=fx,y+d1∇2x,∂ydt=gx,y+d2∇2y.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. The following steps generate the random network and the adjacency matrix element:Step 1: suppose the network consists ofn nodesStep 2: the value of the elementAij of the adjacency matrix is generated as follows: there is an edge between i and j when the random number <p, which is Aij=Aji=1, otherwise Aij=Aji=0. For example, Figure 2 shows a random network structure with a connection probability of p=0.05.Figure 2
Network structure is displayed when connection probabilityp=0.05.Then, the balance equation (3) corresponding to each network node can be rewritten as follows:(11)dxidt=fxi,yi+d1∑jLijxi,dyidt=gxi,yi+d2∑jLijyi,where Lij=Aij−kiδij is the Laplacian matrix, ki is the degree of note i, δij=0 if i≠j, and δij=1 otherwise.The general solution of equation (3) can be expressed as follows:(12)xi=∑m=1Ncmβmeλmtϕim,yi=∑m=1Ncmeλmtϕim,where ∑jLijϕjm=Λmϕim.Firstly, we consider the system’s unique positive equilibrium stateMx¯,y¯ after joining the network diffusion, where fx¯,y¯=gx¯,y¯=0. Suppose that the eigenvalues of the Laplacian matrix is Λii=1,2,…,n, which is arranged in descending order, i.e., 0=Λ1>Λ2>⋯>Λn. Then the Jacobian matrix of each node in the network diffusion system is expressed as follows:(13)Ji=b−1+d1Λia2−b−a2+d2Λi.Then the characteristic function of system (11) can be written as follows:(14)λ2−b−1−a2+d1Λi+d2Λiλ+d1d2Λi2+b−1d2−a2d1Λi+a2=0.From Section2.1, we can know that b−1−a2+d1Λi+d2Λi<0, so the necessary and sufficient condition for Turing instability is that there exists at least one Λi to satisfy wΛi=d1d2Λi2+b−1d2−a2d1Λi+a2<0 or Reλ>0. 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 satisfies 4d1d2a2−b−1d2−a2d12/4d1d2=0 (Figure 3). It is obvious that this equation has two roots k1c2,k2c2, when Reλk2=0. When the continuous system is unstable, Reλk2=0 is established in Figure 4, that is, the value range is k2∈D=λk1c2<λ<k2c2. In addition, the relationship between the network node degree ki and the local eigenvectors −Λi of the Laplacian matrix and the eigenvalues can be obtained in Figure 5.Figure 3
A region diagram of the occurrence of Turing instability for diffusion coefficient.Figure 4
Linear stability analysis. The relationship graph betweenReλ about Λik2 when a=1, b=1.5, p=0.1, and d1=0.01. The critical values are Λc=−2+6/10+46 and d2c=0.1+0.046.Figure 5
The role ofki. (a) The relationship graph between degree ki and eigenvalue −Λi. (b) About the dispersion relationship of ki is shown.
(a)(b)Lemma 2.
For the network Laplacian matrixL, generally kmax=maxki,kmin=minki and Λ is the eigenvalue of L. Therefore, Λ∈C=Λ−2kmax<Λ<M0−kmin.Lemma 3.
For a network-organized system, the system is always stable when all the eigenvalues of the network Laplacian matrix are not in the instability regionΛ∩D=Φ (empty set), the instability occurs when Λ∩D≠Φ, and D∩C≠Φ.The detailed proofs of Lemmas2 and 3 are given in reference [22].


## 2.3. Application of Mean-Field Approximation

Then, 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:(15)dxidt=fxi,yi+d1Cx−kixi,dyidt=gxi,yi+d2Cy−kiyi.where Cx=∑j=1NAijxj,Cy=∑j=1NAijyj, and ki 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. Therefore, when the role of other molecules is fixed, and only the equilibrium stateMx¯,y¯ of the system is considered, the corresponding a single node system can be expressed as follows:(16)dxidt=fxi,yi+d1x¯−kixi,dyidt=gxi,yi+d2y¯−kiyi.From the linear analysis of system (16), the characteristic equation is given as follows:(17)λ2−b−1−a2−d1ki−d2kiλ+d1d2ki2+a2d1−b−1d2ki+a2=0.From the above analysis of the stability of Brusselator, system (16) without the diffusion effect is stable when 0<b<1+a2 holds. Assuming that λ1andλ2 are the eigenvalues of equation (17), the stability of system (16) depends on the sign of wki, where wki=d1d2ki2+a2d1−b−1d2ki+a2. If wki>0, system (16) is stable, otherwise unstable.


## 3. Results

In this section, we present an explanation of theoretical results based on chemical mechanisms. Firstly, the network adjacency matrixA is a symmetric matrix which is randomly generated based on a random network with probability p. The relationship between the Laplacian matrix and the adjacency matrix of the network is Lij=Aij−kiδij, where ki is the degree of the node i and δij is a Dirac function, if i=j, δij=1, otherwise δij=0.We choose the parametersa=1andb=1.5 and derive the relationship between the diffusion coefficients d1 and d2 when the Hopf bifurcation occurs from Figure 3, which shows the critical condition for d1 and d2 when Hopf instability occurs. When diffusion coefficients d1 and d2 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 d2>10+46d1. When the ratio of d2 to d1 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.The bifurcation diagram of the system concerningd2, 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.Figure 6
The bifurcation aboutd2 illustrates the system is stable when p=0.0004 and p=0.69.Firstly, we obtain the bifurcation diagram (Figures7–11(a)) of the reaction-diffusion system with respect to d2 when d1=0.01 under different parameter p, which shows that the bifurcation point is consistent with the critical value for theoretical analysis.Figure 7
(a) The bifurcation diagram aboutd2 when p=0.0005 illustrates that the Turing instability occurs. (b) The pattern formation of the system when p=0.0005 and d2=1.5 illustrates the system is unstable.
(a)(b)Figure 8
(a) The bifurcation diagram ofd2 when p=0.1. The pattern formation of the system when p=0.1. (b) The system is stable when d2=0.1. (c) Turing instability occurs when d2=0.3. (d) When some −Λi fall into the instability region, Turing instability occurs in the system.
(a)(b)(c)(d)Figure 9
(a) The bifurcation diagram ofd2 when p=0.2. The pattern formation of the system when p=0.2. (b) The system is stable when d2=0.1. (c) Turing instability occurs when d2=0.3. (d) When some −Λi fall into the instability region, Turing instability occurs in the system.
(a)(b)(c)(d)Figure 10
(a) The bifurcation diagram ofd2 when p=0.33. The pattern formation of the system when p=0.33. (b) The system is stable when d2=0.1. (c) Turing instability occurs when d2=0.3. (d) When some −Λi fall into the instability region, Turing instability occurs in the system.
(a)(b)(c)(d)Figure 11
(a) The bifurcation diagram ofd2 when p=0.43,. The pattern formation of the system when p=0.43. The system is stable when d2=0.1 (b) and when d2=0.3 (c). (d) When some −Λi fall into the instability region, Turing instability occurs in the system.
(a)(b)(c)(d)In Figure7, we can see that when p=0.0005>0.0004, d2 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 d2. Figure 7(b) shows that the system will appear unstable when d2=1.5 and p=0.0005.Figure8 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 d2=0.1, Figure 8(b) demonstrates that the system is stable. But when p=0.1 and d2=0.3, the corresponding pattern formation (Figure 8(c)) is unstable. When p=0.1 and d2=0.3, the real part of the eigenvalues of the system characteristic equation (17) changes with ki and −Λi is shown in Figure 8(d), which shows that the system is unstable and Turing instability occurs.Whenp=0.2, the bifurcation diagram of xi on d2 is shown in Figure 9(a), which illustrates that the system’s equilibrium point is stable when d2 is lower than the critical bifurcation value. Still, Turing instability occurs when d2 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 d2=0.1 and the system is unstable when d2=0.3, which is consistent with the bifurcation diagram on d2. From Figure 9(d), it can be concluded that when d2=0.3, there exist some ReΛi falling into the unstable region, so the system is unstable.Whenp=0.33, the bifurcation diagram of xi on d2 is shown in Figure 10(a), which illustrates that the system’s equilibrium point is stable when d2 is lower than the critical bifurcation value. Still, Turing instability occurs when d2 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 d2=0.1 and the system is unstable when d2=0.3, which is consistent with the bifurcation diagram on d2. Figure 10(d) shows that when d2=0.3, there are some ReΛi falling into the unstable region, so the system is unstable.Forp=0.43, the bifurcation diagram of concentration xi with respect to d2 can be drawn as Figure 11, and the bifurcation point is consistent with the Turing instability threshold. The bifurcation diagram shows that when d2 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 d2=0.1, but unstable when d2=0.3. In Figure 11(d), there exist some −Λi belonging to D, so Turing instability occurs when d2=0.3.Whenp=0.53, we can give the bifurcation diagram of d2, as shown in Figure 12(a), and the pattern formation when d2=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 d2=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 d2=0.3.Figure 12
(a) The bifurcation diagram ofd2 when p=0.53. The pattern formation of the system when p=0.53. (b) The system is stable when d2=0.1. (c) Turing instability occurs when d2=0.3. (d) When some −Λi fall into the instability region, Turing instability occurs in the system.
(a)(b)(c)(d)But whenp=0.8, the bifurcation of d2 and the corresponding pattern formation of the system are exhibited in Figure 13, demonstrating the system is stable. And the pattern formation, shown in Figures 13(b) and 13(c), indicates that when d2=0.1 and d2=0.3, the system is stable. Figure 13(d) demonstrates that the system’s eigenvalues do not fall in the Turing instability region Λ∩D=Φ, so when d2=0.3, the system is stable. As the connection probability increases, the system’s bifurcation threshold will be increased or even disappeared.Figure 13
(a) The bifurcation diagram ofd2 when p=0.8. The pattern formation of the system when p=0.8. The system is stable when d2=0.1 (b) and d2=0.3 (c). (d) When all −Λi do not fall into the instability region, the system is stable.
(a)(b)(c)(d)From Figure4, we can conclude that Turing instability occurs when ki∈8,39 by the mean-field theory. From Figure 14(a), we can obtain that when ki=7<8,39 does not belong to the instability region, the system will remain in a stable state. When ki=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 ki=50>8,39.Figure 14
The mean-field approximation. (a) The system is stable whenki=7 is not in instability region. The system is unstable when ki=20 (b) and ki=30 (c) are in instability region. (d) The system is stable when ki=50 crosses instability region.
(a)(b)(c)(d)In Figure15, we give the bifurcation diagram of the deterministic system (16) with respect to the node degree. The figure shows that the system has a transcritical bifurcation about the parameter ki near the equilibrium point, where the bifurcation points are, respectively, k1c=7.033 and k2c=37.86, further confirming the system instability region about ki.Figure 15
The bifurcation diagram ofki. The black line indicates that the equilibrium point of the system is unstable, and the red line indicates that the equilibrium point is stable.For how the connection probabilityp 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. The approximate instability region of the reaction-diffusion system for the Brusselator model with respect to d2 and p is shown in Figure 16.Figure 16
The instability region is related to the parametersd2 and p, where the red region is the instability region.Therefore, from Lemma3 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 ofd2 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/n2<p<d2c/d14/5lnn/n, where 1/n2 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.


## 4. Conclusion

In conclusion, firstly, the Turing instability critical value of the Brusselator model with a random network is obtained in Figure3. 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/n2<p<d2c/d14/5lnn/n. When the control parameters and d1 are given, the bifurcation of the reactant concentration on d2 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. The 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. Therefore, 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 reaction-diffusion system with respect to the network node degree is obtained. Whenki=7<8,39, the system remains stable at the equilibrium point. When ki=2030∈8,39, the system has Turing instability. When ki=50>8,39, the system is stable at the positive equilibrium point. And the bifurcation of ki 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.



---

*Source: 1572743-2020-10-05.xml*