Spatiotemporal patterns in a general networked activator–substrate model

To understand the spatiotemporal pattern formation in the random networked system, a general activator–substrate model with network structure is introduced. Firstly, we investigate the boundedness of the non-constant steady state of the elliptic system of the continuous medium system. It is found that the non-constant steady state admits their upper and lower bounds under certain conditions. Then, we establish some properties and nonexistence of the non-constant steady state with the no-flux boundary conditions. The main results show that the diffusion rate of the activator should be less than that of the substrate. Otherwise, there might be no pattern formation of the system. Afterward, a general random networked activator–substrate model is presented. Conditions of the stability, the Hopf bifurcation, the Turing instability and the co-dimensional-two Turing–Hopf bifurcation are yielded by the method of stability analysis and bifurcation theory. Finally, a suitable sub-system of the general activator–substrate model is chosen to verify the theoretical results, and full numerical simulations have well verified these results. Especially, an interesting finding is that the stability of the positive equilibrium will switch from unstable to stable one with the change of the connection probability of the nodes, which is different from the pattern formation in the continuous medium systems.


Introduction
In the real world, reaction-diffusion is one of the most basic movement processes, and can be used to describe the spread of diseases, chemical reactions, biological migration and other phenomena. Hence, many coupled reaction-diffusion equations are used to understand this process in a continuous space. An interesting and challenging research area is to understand the Turing patterns by a coupled reaction-diffusion equation [1][2][3][4]. Along this way, there have been plenty of accumulated achievements to enrich pattern formations within coupled reaction-diffusion systems, for example, the stripe patterns, the stripe and the spotted mixed patterns in a one-dimensional region [5,6], and more complex hexagonal patterns, the labyrinthine-like stripe patterns, Eckhaus patterns, chaos in a two-dimensional region [7,8]. What is more, the spatiotemporal patterns also have been investigated where the Turing mode and the non-Turing one intersect in reaction-diffusion equations. For instance, a codimension-two Turing-Hopf bifurcation was considered in [9,10] by the technique of the normal form theory and the multiple time scales, respectively.
It is noticed that models in the references mentioned above are governed by the continuous systems. However, as it is known that in many cases the reactiondiffusion process under scrutiny is defined in a discrete medium, rather than in a continuous one, such as disease spread and control, population dynamics, and so on. Therefore, could such rich spatial and spatiotemporal behaviors happen in the discrete reaction-diffusion systems? Keep all in mind, a reaction-diffusion model with the network structure could be well worth considering. Asllani et al. [11] showed that the homogeneous fixed point can become unstable due to the topology of the network, which cannot be induced on undirected graphs, and the traveling waves or quasi-stationary patterns could be found in a network system. Diego et al. [12] provided a first general theory of Turing network topology, proving why three key features of a Turing system are directly determined by the topology. They also yielded that in experimental systems, it was easier to obtain reliable information about the topology of a network than precise values of diffusion constants. The instability and pattern formations on several different networks including deterministic networks and random networks with time delay were explored by Chang et al. [13], and wave patterns could be generated on different networks by numerical simulations. The Hopf bifurcation and its detailed nature were considered in a network-organized semiarid vegetation model with time delay by Tian et al. in [14]. For more dynamical results on the network-organized systems, one could refer to Refs. [15][16][17][18] and references cited therein.
However, note that the results about the influence of the connection probability of nodes are few [19]. Inspired by this discovery, we will explore the influence of the connection probability of nodes in a networked model. To this end, we first introduce a continuous medium system of the form where u = u(x, t) and v = v(x, t) are the concentrations of the activator and the substrate, respectively; constants d 1 and d 2 stand for the diffusion coefficients of activator u and substrate v, respectively; is recorded as the Laplacian operator; k ≥ 2 and f (v) is a monotone increasing function with respect to v; spatial domain is a bounded domain in R n (n ≥ 1) with the boundary ∂ ; n is the outward unit normal vector along ∂ and ∂n represents the operator of directional derivative along the direction n; parameters c, β, d 1 and d 2 are positive constants from the point of view of chemistry. The dimensionless system (1) is commonly called the depletion model or the activator-substrate model, which was firstly proposed by Gierer and Meinhardt in 1972 [20]. Meinhardt and Klingler [21] employed this model to explore pigmentation and relief-like patterns on mollusk shells when f (v) = v and replaced u k by u 2 1+κu 2 +ρ 0 , where κ > 0 and ρ 0 represents a small basic production of the activator. A similar result was carried out by Fowler et al. [22] to describe the pigmentation patterns on shell surfaces. The Hopf bifurcation and its direction, as well as the Turing instability, were studied by Wu et al. [23], with the monotone increasing function f (v) = v and k = 2. Also, the asymptotic stability of constant steady states and the steady state bifurcations from constant steady states were investigated in both one-dimensional and two-dimensional kernel cases by Wang et al. [24]. For more results about the activator-substrate model (1), one can refer to Refs. [25,28]. Now, assume that each node of the network owns the homogeneous density in a small community, then a general networked activator-substrate model takes the form Here, we assume that the networked model (2) is defined on an undirected network with N nodes and there are no self-loops; u i and v i are the concentrations of the activator and the substrate on node i, respectively; is the N × N Laplacian matrix of network, set to be i j = k i δ i j − L i j , where L is the adjacency matrix encoding the network connection, namely it satisfies L i j = 1 if there is a link connecting from patch i to patch j, otherwise, L i j = 0 if there is no any link connecting from patch i to patch j, and δ i j is the Kronecker's delta. In addition, the degree of the ith node is defined by k i = N j=1 L i j , and p(0 ≤ p ≤ 1) is set to be the connection probability between node i and node j for i = j.
In this work, we shall investigate the dynamics of the general activator-substrate model (1) and networked model (2), respectively. More precisely, the boundedness of the positive non-constant steady state will be studied for the elliptic system of the continuous general activator-substrate model (1). With the help of the maximum principle and Harnack's inequality, one shows that it admits the upper bounds and the lower bounds for any positive solution (u(x), v(x)) under certain conditions. Also, we will explore the sufficient conditions to ensure the nonexistence of the non-constant steady state utilizing the Poincaré's inequality and Cauchy-Schwarz inequality. It is found that there is no nonconstant steady state of the elliptic system if the diffusion rate of the activator is greater than the diffusion rate of the substrate. This is very useful information for studying the formation of spatiotemporal patterns. Moreover, we would like to mention that all the results obtained above employ the homogeneous Neumann boundary conditions and assume that f (v) satisfies (C 0 ), where For the networked activator-substrate model (2), some general results about the existence of internal positive equilibria and their stability, the Hopf bifur-cation and the Turing-Hopf bifurcation are given, and a numerical example is considered to verify these general results. As we know, the Turing-Hopf bifurcation is a codimension-two bifurcation, where both the Hopf mode and the Turing mode will intersect. Hence, the formation of patterns may be more complex than that of patterns near the Hopf mode and the Turing mode. Compared with the works done in [14] and references cited therein, this work has two highlights: (i) In addition to the Hopf bifurcation and the Turing instability, the spatiotemporal patterns near the Turing-Hopf bifurcation are also considered theoretically and numerically. (ii) The numerical simulation results indicate that the connection probability p has a crucial influence on the Turing pattern formation, namely the stability of the internal positive equilibrium will switch from unstable to stable one with the changing of the connection probability p, although the parameters of the model (2) are fixed in the Turing instability region.
The layout of this paper is structured as follows. In Sect. 2, the boundedness and the nonexistence of the non-constant steady state of the related elliptic system (1) are carried out. In Sect. 3, stability, the existence of the Hopf bifurcation, the Turing instability and the Turing-Hopf bifurcation of the general networked activator-substrate model (2) are considered. Then, a sub-system of the general networked activatorsubstrate model (2) is introduced in Sect. 4, and full numerical simulations are carried out in Sect. 5 to verify the theoretical validity. Finally, some discussions and conclusions are made in Sect. 6.

Non-constant steady state
In this section, we shall investigate the nonexistence of the non-constant steady state of the elliptic system; it is governed by

A priori estimates
and ω( and where positive constant c * only depends on Proof We assume that (u(x), v(x)) is a positive solution of system (3), and By Lemma 1 and the second equation of (3), we have As a result, by the first equation of (3), we have On the other hand, define a linear function of the form then it follows that In addition, we assume that ω(x) has maximum at x = z 1 , i.e., ω(z 1 ) = max x∈¯ ω(x). Then by means of Lemma 1 and Eqs. (6) and (7), one obtains Consequently, we have It is noticed that v(x) ≤ c, then one has Therefore, by Lemma 2 we claim that there is a constant By employing Lemma 1 and the second equation of (3) again, we have The proof is completed.
In what follows, we investigate nonexistence of the non-constant steady state solution of system (3). To do so, let 0 = λ 0 < λ 1 < λ 2 < · · · < λ i < · · · and lim i→∞ λ i = ∞ be the complete set of eigenvalues of the operator − with no-flux boundary conditions in .

Nonexistence of non-constant steady state
Suppose that (u(x), v(x)) is a positive solution of system (3), we denote their averages over domain bȳ Then, it is not difficult to verify a fact that Let φ = u −ū and ϕ = v −v, then we know that φdx = 0 and ϕdx = 0.
is an arbitrary positive solution of system (3), then by Cauchy-Schwarz inequality and integrate by parts, one obtains By employing the Poincaré's inequality, that is where λ 1 is the first positive eigenvalue of − over domain with respect to the no-flux boundary conditions, we therefore obtain This leads to By employing the Poincaré's inequality again, one achieves This completes the proof. Now, we consider nonexistence results of the nonconstant steady state solution of system (3).

Lemma 5
Assume that c, β, d 1 and d 2 are positive constants in (3), then Now we yield the following result.
) be a non-negative steady state solution of system (3). Then, multiplying the first equation of (3) by φ and integrating by parts over domain , one yields where we denote By Lemma 5, one has Hence, the Poincaré's inequality and Cauchy-Schwarz inequality reveal that Similarly, multiplying the second equation of (3) by ϕ and integrating by parts over domain , we have As a result, one deduces where d * 1 and d * 2 are set to be holds, then one has ∇φ = ∇ϕ = 0. This indicates that (u(x), v(x)) must be a positive constant solution of elliptic system (3).

constants in (3) and ⊂ R n is a bounded domain with sufficient smooth boundary, then elliptic system (3) has no non-constant steady state solution if
Proof multiplying Eq. (7) by φ and integrating it by parts, we have From Eq. (6), one obtains We therefore obtain Combining Eqs. (8) and (9), we have Then, by employing Eqs. (6) and (10), one has Therefore, the Poincaré's inequality and Cauchy-Schwarz show that It then leads to Consequently, we claim that (u(x), v(x)) must be a constant solution of elliptic system (3) if The desired result is confirmed.
Remark 1 By means of Theorem 3, we find a fact that This shows that system (3) always has no non-constant solution if d 1 is large enough and d 2 is fixed. Hence, given the spatial or spatiotemporal patterns of the system (3) we would restrict d 1 < d 2 , namely the diffusion rate of the activator should be less than that of the substrate. Otherwise, there might be no pattern formation of system (3).

Spatiotemporal patterns in the networked model
In this section, we shall investigate the Hopf bifurcation, Turing instability and the Turing-Hopf bifurcation of the networked system (2). It is easy to check out that system (2) has a positive equilibrium E * = (u * , v * ) = (u * , c − βu * ) when 0 < u * < c β . Keeping this in mind, we now perform the linear stability analysis of system (2) near internal positive equilibria E * . To do so, let δu i = u i − u * and δv i = v i − v * , then the linear system of (2) evaluated at E * can be described as follows where Let 0 = 1 > 2 > · · · > N be the eigenvalues of the discrete Laplacian matrix , and we assume that is the subspace generated by the eigenfunctions associated to the eigenvalue i . As a result, the general solution of system (2) can be written as follows with N j=1 i j φ j = i φ i . Putting (12) into (11), then for each node i = 1, 2, · · ·N , we have It follows that Hence, we have the characteristic equation for each node with the form where . Then, the characteristic roots of Eq. (13) can be exhibited of the form Now, we have stability results of the internal positive equilibrium E * as follows.

then networked system (2) undergoes a Hopf bifurcation, where
Proof Let's recall Eqs. (13) and (15). The stability of the positive equilibrium E * can be determined by the sign of T i (d 2 , β) and D i (d 2 , β). That is, the positive equilibrium E * is locally asymptotically sta- This implies there is at least one eigenvalue of the characteristic Eq. (13) with positive real part, so E * is unstable, no matter what the sign of T 0 (d 2 , β) is. Conclusion (i) is valid. Condition 0 < k k−1 < β H indicates D 0 (d 2 , β) > 0. So it is readily to verify that (ii) and (iii) are true. On the other hand, taking the derivatives with respect to β yields Therefore, the Poincaré-Andronov-Hopf bifurcation theorem shows that the local system of (2) undergoes the Hopf bifurcation at E * when β = β H . This completes the proof.
In what follows, we focus on the Turing instability of the positive equilibrium E * when diffusion is presented.
Theorem 5 Suppose that c, β, k ≥ 2 are positive con- Proof From Theorem 4, we know that E * is locally asymptotically stable when 0 < k k−1 < β H and β < β H . This leads to T 0 (d 2 , β) < 0 and D 0 (d 2 , β) > 0. On the other hand, note that 0 It is noticed that F u > 0 and G v < 0, we thus obtain So if then the network system (2) becomes Let α = c β . Then, when α = 2, the system (18) has a unique internal positive equilibrium E 1 = (1, β). When α > 2, the system (18) has two internal positive equilibria E 2 and E 3 , where From [23], the internal positive equilibrium E 1 is non-hyperbolic and E 2 is an unstable saddle. Therefore, we shall mainly focus on the dynamical behaviors of the internal positive equilibrium E 3 , denoting it by E * , that is, Obviously, we have 0 < u * < α. Then, we obtain the fact that As a result, we exhibit the main results of system (18)  Proof It is easy to obtain Then based on Theorem 4, the results are valid.
Proof Similar to the proof of Theorem 5, here we omit the detailed proof process.
Proof By the same manner in Theorem 6, one could show these results. Here, we omit the detailed proof process.

Numerical simulations
In Sect. 4, we perform the stability analysis and bifurcation analysis of the networked system (18) and the related results have been exhibited in Theorems 7-9. Now we shall carry out some numerical simulations to verify these results. We fix N = 100, that is, the numbers of the nodes in the random network system (2) are 100. Moreover, we set p(0 ≤ p ≤ 1) to be the connection probability of different nodes u i and u j . Firstly, to exhibit the distribution of the solutions u i and v i of the system (2), we assume that the network is a Watts-Strogatz type [29]. From the point view of mathematics, a Watts-Strogatz network is an undirected graph. Some of the connections between the nodes are determined and some of them are random. It possesses the homogeneous degree distribution [14]. It is noticed that the distribution of solutions u i and v i under Watts-Strogatz network is very similar, so we only perform the distribution of u i under Watts-Strogatz network. Here, parameters in random network system (2) are chosen as c = 2.2, β = 1, d 1 = 1.5 and d 2 = 10, the average degree is 4 and the connection probability p is a variable. It is found that u i and u j will eventually become a closed circle when connection probability p = 0, where 1 ≤ i, j ≤ 100 and i = j. However, when one chooses 0 < p < 1, such closed circle will be destroyed. More precisely, u i will randomly walk into different positions at the same time t. Hence, connection probability p has a vital influence on the distribution of u i , see Fig. 1 for more details. Now we choose c = 23.58, β = 7.86 and connection probability p = 0.05 then one has α = 3, β H =   then we have Re(λ i ) < 0 and D i (d 2 , β) > 0. That means that no Turing pattern will emerge in the random network system (2). Figure 4 shows the internal positive equilibrium E * is locally asymptotically stable. Here, we choose d 2 < d T 2 and the other parameters are fixed in Fig. 3. However, when taking diffusion coefficient d 2 = 10 > d T 2 , we find that the Turing instability emerges and such instability induces the stripe patterns, see Fig. 5. Therefore, the theoretical analysis in Theorem 8 is well verified by the numerical simulation results exhibited in Figs. 4 and 5. Figure 6 displays the bifurcation diagram in the plane of β − d 2 . It is found that the plane of β − d 2 is divided into four regions, namely the stability region, the Hopf instability region, the Turing instability region and the Turing-Hopf bifurcation region, respectively. Obviously, the parameters in Figs. 4 and 5 are located in the stable region and the Turing instability region of Fig. 6, respectively. Now we choose c = 23.58, d 1 = 2 and p = 0.02, then we obtain the fact that d 2 = d * 2 = 10.56 and β = β * = 7.85 with critical value i = ic = −1. This implies that the random networked system (2) may To do so, we take parameters in the Turing-Hopf region, i.e., d 2 = 10.6 and β = 7.86. As a result, the random networked system (2) admits the Turing-Hopf bifurcation and the spatially inhomogeneous periodic solution can be induced by such codimension-two bifurcation, see Fig. 7.
In short, all the theoretical results in Theorems 7-9 are valid.
To further understand the role of connection probability p in the Turing pattern formation. In what follows, we shall fix the parameters in Fig. 5, that is, c = 2.3, β = 1, d 1 = 0.5, d 2 = 10, but change the connection probability p. We only perform the patterns of u i for convenience.
In Fig. 8, we take the connection probability p = 0. The numerical simulation result shows that there are no eigenvalues − i of the discrete matrix located in the Turing instability region. So the internal positive equilibrium E * is asymptotically stable, namely the random networked system (2) will not exhibit the Turing instability.
In Fig. 9, we shall choose the connection probability p = 0.02. The numerical simulation shows that there is some − i in the Turing instability region, then the real parts of characteristic roots are positive, i.e., Re (λ i ) > 0. Therefore, the random networked system (2) undergoes the Turing instability.
In Fig. 10, we shall choose the connection probability p = 0.05. The numerical simulation shows that there is some − i in the Turing instability region. However, compared with Fig. 9, there are fewer − i in fewer and the Turing instability will eventually happen in the random networked system (2). In Fig. 12, we take the connection probability p = 0.09, then we find that there is no − i lying in the Turing instability region. That is, although the same parameters are fixed in Fig. 5, the Turing pattern will not appear in the random network system (2). The right picture of Fig. 12 well verifies such prediction. It is found that the internal positive equilibrium E * becomes a stable one again.
Summarizing the numerical results exhibited in Figs. 8, 9, 10, 11 and 12, we note that the connection probability p has an essential influence on the pattern formation in the random networked system. This is very different from the pattern formation in the continuous medium systems. Precisely, for the continu- There is no eigenvalues of the discrete Laplacian matrix located in Turing instability region, and positive equilibrium E * becomes stable again. Here, c = 2.3, β = 1, d 1 = 0.5, d 2 = 10 and p = 0.09 ous medium system, if the parameters are fixed in the Turing instability region, then the eigenvalues of the Laplacian matrix are fixed. Hence, the Turing pattern may exist. However, for the random network system, although the parameters are fixed in the Turing instability region, the eigenvalues of the discrete Laplacian matrix, i , are not fixed, namely they may be in or out of the Turing instability area with the change of the connection probability p. Hence, apart from diffusion coefficients, the connection probability p also has a crucial influence on the pattern formation of random network systems compared with the continuous medium systems.

Conclusions
In this work, a general activator-substrate model with network structure is introduced to understand the formation of spatiotemporal patterns. We first consider the general continuous activator-substrate model. Here, the maximum principle and the Harnack's inequality are used to investigate the boundedness of the nonconstant steady state of the elliptic system. It is found that the non-constant steady state admits their upper and lower bounds with certain conditions, see Theorem 1. Then, the Poincaré's inequality, Cauchy-Schwarz inequality and the technique of integrating by parts are employed to study some properties and nonexistence of the non-constant steady state with the no-flux boundary conditions. One shows that the general activatorsubstrate model could admit the non-constant steady state if the diffusion rate of the activator is less than that of the substrate, see Theorems 2 and 3. Afterward, a general random network activator-substrate model is introduced. The existence conditions of the stability, the Hopf bifurcation, the Turing instability and the Turing-Hopf bifurcation are yielded by the method of linear analysis and bifurcation theory, see Theorems 4-6, respectively. Finally, a suitable sub-system of the general activator-substrate model is investigated and numerical simulations have well verify these results. Especially, the unstable positive equilibrium E * will become stable with the change of the connection probability p by comparing Figs. 8, 9, 10, 11 and 12. Hence, the connection probability p of the nodes is also an important factor for the formation of the Turing pattern compared with the continuous medium systems. In addition, whether the degree of nodes affects the pattern formation in the random network systems is also worth exploring. That will be the future work.