Cross-diffusion induced spatiotemporal patterns in Schnakenberg reaction–diffusion model

Schnakenberg system is a typical mathematical model to describe activator-depleted kinetics. In this paper, by introducing linear cross-diffusion into Schnakenberg system, we derive cross-diffusion-driven Turing instability conditions. It has been revealed that it is no longer necessary to have long-range inhibition and short-range activation for Turing instability with the help of cross-diffusion. Then, the multiple scales method is applied to obtain the amplitude equations at the critical value of Turing bifurcation, which helps us to derive parameter space more specific where certain patterns such as hexagon-like pattern, stripe-like pattern and the coexistence pattern will emerge. Furthermore, the numerical simulations in both Turing instability region and Turing–Hopf region provide an indication of the wealth of patterns that the system can exhibit. Besides, different initial conditions are employed to help better understanding the complex patterns.


Introduction
Pattern formation can be induced by diffusions which has been proposed by Turing [1] in 1950s. These kinds of pattern therefore can also be called Turing pattern. The two crucial conditions for Turing pattern are: two coupled mechanisms which exist in the reaction system: autocatalytic kinetics and diffusion; the diffusion rate of one morphogen is much larger than the other. Reaction-diffusion systems have been generally investigated in the field of ecology, demography and also chemistry. In 1979, Schnakenberg proposed a chemical reaction model as an example to exhibit limit cycle behavior [2]. The non-dimensional form of Schnakenberg model as a typical activator-depleted substrate chemical reaction-diffusion model is governed by Here, u(x, t) and v(x, t) are the concentrations of two interacting morphogens and d 1 , d 2 are the corresponding diffusion rates; p and q are positive constants representing the chemical sources.
Cross-diffusion first was proposed by Kerner in view of the possible bias of the motion of one morphogen toward or away from another [25]. Later, a competitive population system was proposed by Shigesada et al. involving the idea of cross-diffusion [26]. Since the pioneer works, the concerns on cross-diffusion have attracted the attention of scholars in many fields. Crossdiffusion has been demonstrated to be quite significant for such reasons: cross-diffusion can induce pattern transitions, even in relatively simple systems; crossdiffusion can promote the density of morphogens; cross-diffusion can trigger spatial patterns which is however impossible under the same conditions only with self-diffusion and so on. Recently, cross-diffusion driven Turing pattern, as an interesting and challenging problem, has been subjected to intensive studies [27][28][29][30][31][32][33][34][35][36].
Considering cross-diffusion on Schnakenberg model, the most relative papers are as follows. The conditions for occurrence of Hopf, Turing and wave instabilities of Schnakenberg model with cross-diffusion have been performed through linear stability analysis, and one-dimensional numerical simulations have been illustrated on the basis of weakly nonlinear analysis [37]. Madzvamuse et al. [38] have numerically solved the Schnakenberg model using finite element methods and illustrated how cross-diffusion induces the formation of pattern.
In [39], several processes which can give rise to cross-diffusion have been described such as chemotaxis, ecology, epidemiology and physicochemical reaction-diffusion systems. For chemotaxis, it is necessary to introduce cross-diffusion because that cells direct their motion toward or away from higher concentrations of chemoattractants or chemorepellants. Based on the previous works of Schnakenberg model with cross-diffusion [37,38], we now concentrate on the following Schnakenberg model with cross-diffusion: where d 11 and d 22 are referred as self-diffusion coefficients, d 12 and d 21 are called cross-diffusion coeffi-cients. The constants p and q still represent the chemical exogenous sources.
The main object of this paper is to investigate the effect of cross-diffusion on pattern formation in Turing instability region and Turing-Hopf region when added into nonzero initial conditions: and zero-flux boundary conditions where L denotes the size of the spatial variables x and y, n is the outward unit normal vector of the boundary ∂ . We expect to make the following contributions to Schnakenberg model: • We demonstrate that it is no longer necessary to have long-range inhibition and short-range activation as premise for pattern formation. That is shortrange inhibition and long-range activation can also motivate Turing instability. • For the same parameter values, we take different initial conditions to contrast the obtained patterns. Then, a conclusion has been arrived that under these different initial conditions, the patterns' constitutions keep pace with each other although the distributions of spot and stripe patterns somewhat different. • We attempt to make numerical simulations to obtain patterns in Turing-Hopf region and have observed different spatial patterns.
The rest of the paper is as follows. In Sect. 2, by linear analysis, we obtain the conditions for stability and Turing instability of the positive equilibrium. In Sect. 3, we investigate pattern transition by deriving amplitude equations using multiple-scale analysis. In Sect. 4, some numerical simulations are shown to illustrate the important impact of cross-diffusion on pattern formation. Finally, some conclusions are drawn.

Turing instability analysis
In this section, we consider the influence of crossdiffusion on the stability of the positive equilibrium and investigate cross-diffusion-induced Turing instability of system (2). By simple calculations, the unique coexistence equilibrium of system (2) The Jacobian matrix corresponding to (u * , v * ) is as follows: We seek solutions in the form of with k · k = k 2 , k and λ are the wave number and frequency, respectively. i is the imaginary unit satisfying i 2 = −1, and r = (u, v) is the spatial vector in two-dimensional space and c.c. stands for the complex conjugate terms. Substituting it into the linearization of (2) leads to the following dispersion relation, which gives the function of eigenvalue λ at the wave number k: where The roots of Eq. (4) can be solved as: We choose q as the bifurcating parameter to investigate the occurrence of Hopf bifurcation and Turing instability. Hopf bifurcation occurs when where λ k is the root of the characteristic equation (4). Namely Hopf bifurcation curve is defined as in the p − q parameter plane from which we can obtain the critical value q H for Hopf bifurcation. For Turing instability to be achieved and different spatial patterns to form, we should ensure at least one root of the characteristic equation (4) with positive real part in the parameter region where system (2) without diffusion terms at E * is asymptotically stable. Re(λ k ) will be positive for certain k = 0 at which J k < 0. That is to say that Turing bifurcation occurs when where k T stands for the critical wave number of Turing bifurcation. The steady state is marginally stable In p − q plane, let us discuss the bifurcations represented by the formulas (5) and (6). Figure 1 shows two cases for bifurcation graph. In Fig. 1a, the unstable region consists of Turing instability region and Turing-Hopf region for d 11 Figure 2 shows that the coexistence equilibrium E * is stable corresponding to the non-diffusion system when q > q 1 = q H , and becomes Turing instability when added the diffusive effects when q < q 5 = q T . In the following section, amplitude equations are extended for q ∈ (q 1 , q 5 ) = (0.7789, 1.5226).

Amplitude equations
From the analysis above, system (2) will experience Turing instability at q = q T for the critical wave number k = k T . In order to further divide Turing instability region in p − q plane into certain parts where various Turing patterns will arise, we employ the method of multiple scales [40] to compute the amplitude equations taking q as the bifurcating parameter. In the following, we describe the computation of amplitude equations via three steps.
Step 1 Rewrite out the linearized form of model (2) at the positive equilibrium point (u * , v * ) as follows: satisfying the following equalities respectively, among which L + T is the adjoint operator of L T .
Step 2 Close to the onset q = q T , we expand the solutions of model (2) as the following asymptotic expression: and Step 3 By the above expressions, together with the amplitude equations about A j are as follows: where with

Spatial dynamics
Associated with the decomposition A j = ρ j exp(iθ j ), (7) can be converted to the following equations with the mode ρ j =| A j | and a corresponding phase angle θ j by separating the real and imaginary parts: with θ = θ 1 + θ 2 + θ 3 . We expect different kinds of patterns of (2) corresponding to different constant solutions to (9). The existence and stability of the constant solutions can be concluded as follows.
(IV) The mixed states (MS): the solution is with g 2 > g 1 , which exists when μ > μ 3 and is always unstable.
Taking  (8), we obtain that h = −1.1639, g 1 = 7.4554 and g 2 = 21.0932. Close to q T , then we can get μ 1 = −0.0068, μ 2 = 0, μ 3 = 0.0543 and μ 4 = 0.2622. In Turing instability region, in other words for 0.7789 < q < 1.5226, we further divide Turing region into different parts where various spatial patterns arise. Summarizing the above four circumstances, we show the results intuitively by Fig. 3. For q H < q < q T , which indicates μ 1 < μ 2 < μ < 0.4884, the stationary state of the diffusive system (2) begins to lose stability and first the hexagon patterns through the coexistence equilibrium begin to form. In this parameter region, the sign of h according to Eq. (8) determines whether the hexagon patterns are from H 0 (i.e., θ = 0) or H π (i.e., θ = π ). Considering h < 0 according to (8), the hexagon pattern is of H π . The emergence of the stripe patterns derives from supercritical bifurcation, and they are only stable for μ > μ 3 . Then, q continues to grow, and when μ satisfies μ 1 < μ 2 < μ 3 < μ < μ 4 , there exists a bistable state of the coexistence of hexagon patterns and stripe patterns. The system evolves to totally stripe patterns for μ satisfying μ 1 < μ 2 < μ 3 < μ 4 < μ.

Numerical simulations
For the remainder of this section, we present some numerical simulations on system (2) through the finitedifference method constructed in two spaces dimension via MATLAB by Garvie [41]. This algorithm is stable and convergent with time step t = 0.01. For the twodimensional approximations, a uniform subdivision of the discrete domain = [0, 100] × [0, 100] with the grid points (x i , y i ) = (i h, i h), i = 0, 1, . . . , 200 with the space step h = 0.5 has been adopted.
From an ecological point of view, the densities of populations or morphogens are oscillatory based on the situation of primary interest [41]. Thus, in the numerical simulations, there are two kinds of initial data most commonly used. One is the random spatial distributions, and the other is small spatially homogeneous perturbation of the steady state of E * in view of special need. Through the numerical simulations, the distributions of morphogens u and v observed are always of the similar structure. As a result, we only show the pattern formation about morphogen u for instance. Through the numerical calculation, the coefficients g 1 , g 2 , h, μ 1 , μ 2 , μ 3 , μ 4 in amplitude equations (9) can be obtained. As parameter q varies, the corresponding control parameter μ lies in different bifurcating domain from Fig. 3. Here, we show three kinds of categories of Turing patterns. In every pattern, the yellow (blue) represents the higher (lower) concentration of morphogen u.

Pattern formation for Turing instability region
Before the random disturbances, the system initially is at the stationary state of E * , and we run the simulations until they reach a stationary state or the pattern does not change any more with time elapsing.
The initial conditions for Figs. 4, 5, 6 and 7 are taken as a uniformly distributed random perturbation in the order of 0.0001 around the steady states formally as:  Figure 4 shows the formation of the spatially hexagon-like patterns of morphogen u, respectively, after 0 iteration, 60,000 iterations, 68,000 iterations, 78,000 iterations, 85,000 iterations and 136,000 iterations for q = 1.5. Then, μ = 0.0148 and μ 1 < μ 2 < μ < μ 3 < μ 4 . Initially, the random distribution leads to the irregularly transient pattern in the domain (cf. Fig. 4a). Afterward, the H π hexagon patterns format very slowly, and until about 60,000 iterations, some regular spots begin to emerge (cf. Fig.  4b). As the density of morphogen u changes, the random perturbations lead to the gather of hexagons, and the pattern becomes clearly whether from the color or from shape (cf. Fig. 4c-e). After the evolution, the pattern ends with uniform hexagons like "honeycomb" in dense arrangement and the structure does not change any more until 136,000 iterations (cf. Fig. 4f). The final pattern consists of bluish green hexagons on the yellow background, which we can call them "cold hexagons." From chemical point of view, q = 1.5 means more addition exogenous substances compared with q = 1.4 in Fig. 5 and q = 1.3 in Fig. 6. Along with the increase in v, the growth of u has been lessened. Thus, the diffusion and decay of u dominate over the production of u far from regions of high concentration of u. Figure 5 shows the formation of the stationary mixed hexagon and stripe patterns of morphogen u, respectively, after 0 iteration, 10,000 iterations, 18,000 iterations, 30,000 iterations, 45,000 iterations and 58,500 iterations for q = 1.4. Then, μ = 0.0805 and μ 1 < μ 2 < μ 3 < μ < μ 4 . Initially, the random distribution leads to the irregularly transient pattern in the domain (cf. Fig. 5a). Afterward, the irregularly transient pattern changes fast due to the change of the density of morphogen u and some regular hexagons begin to emerge (cf. Fig. 5b). Some of the hexagons continue to aggregate and obvious stripes spring up which leads the hexagons together with the stripes prior formed gradually to stay apart from each other (cf. Fig. 5c). After 30,000 iterations, a few stripes emerge competing with the remaining part of the hexagons (cf. Fig.  5d, e). After the evolution little by little, the hexagons and stripes patterns are staggered in the whole domain,  Fig. 5f). Figure 6 shows the formation of the stationary stripe patterns of morphogen u, respectively, after 0 iteration, 10,000 iterations, 20,000 iterations, 30,000 iterations, 80,000 iterations and 117,000 iterations for q = 1.3. Then, μ = 0.1462 and μ 1 < μ 2 < μ 3 < μ < μ 4 . Model dynamics exhibit a transition from hexagonstripe growth to stripe replication, i.e., hexagons decay and stripe patterns appear. Starting with the random distribution, the irregularly transient pattern (cf. Fig.  6a) changes fast to the mixture of hexagon-stripe pattern (cf. Fig. 6b). The later random perturbations make these hexagons vanish (cf. Fig. 6c-e), ending with the time-independent stripe pattern in the whole domain (cf. Fig. 6f). In this case, the numerical simulations cannot correspond to the theoretical analysis due to that the control parameter q is far away from the critical Tur-ing bifurcation value q T = 1.5226. What is more, these kind of stripe patterns can be called as "labyrinthine" pattern. According to the study from Cartwright [42], the labyrinthine patterns stem from Turing instability from an initially homogeneous system by symmetry breaking. In the cell cultures, the labyrinthine patterns come from the warfarin treatment [43].
When the control parameter q locates further away from q T = 1.5226 in Turing instability region, we show the snapshots of morphogen u at different values of q with d 11 = 0.1, d 12 = 0.1, d 21 = 0.1, d 22 = 2 and p = 0.1. For q = 1.25, one can see that the stripe-like patterns almost prevail over the whole domain only with several spots (cf. Fig. 7a). When q decreases to 1.2 and 1.05, the spotted patterns have been well-distributed both in shapes and colors in the whole domain. This kind of spot pattern seems different from that of Fig. 4f and the color varies from bluish  Fig. 7b, c). This kind of spot pattern consists of yellow hexagons on the dark blue background, which we can call them "hot spots." Throughout the numerical simulations of morphogen u above, we find that with the decrease in control parameter q, the pattern sequence "cold hexagons → mixture of cold hexagons and stripes → stripes (labyrinthine) → mixture of hot spots and stripes → hot spots" is observed. In a word, the pattern changes from hexagon to stripe and eventually to spot. The appearance of hexagon pattern is derived from that the fluctuation of the concentrations of activators u fires the central peak of activity with a substantial spread of the activators by diffusion and regularly may be limited by the morphogen v. Then, with the diffusion of the activator u, the pattern changes from hexagon to stripe. Because that diffusion can cause the activated morphogens to have the tendency to appear in coherent patches. In stripes, activated morphogens have activated neighbors along the stripe without non-activated morphogens close by, in which the morphogen v can be dumped. Then, the increase in v will lessen the growth of u and spot patterns show up. Due to the limit owing to saturation, the activated morphogens remain in a scattered arrangement with the activation at a lower level. That is the reason why the cold hexagon and hot spot patterns seem unlike.
In Fig. 8, all the parameters are fixed the same as Figs. 5 and 7a. When the initial values change into: and u(x, 0) = u * + 0.001 × (cos(x) + cos(y)), v(x, 0) = v * + 0.001 × (cos(x) + cos(y)) (11)   (10); b the initial condition (11) separately, the pattern's constitutions of mixture of spots and stripes keep pace with Figs. 5 and 7a, while the structures in which spots patterns and stripe patterns coexist are different from those in Figs. 5 and 7a due to the fact that the symmetric initial conditions are employed instead of the initial conditions with uniformly distributed random perturbation.

Pattern formation in Turing-Hopf region
In this subsection, we perform some numerical simulations in Turing-Hopf region. To this end, we shall set parameters as follows: p = 0.13, q = 0.2, and then we obtain some interesting simulation results. In Fig. 9, the control parameters p and q are set as p = 0.15 and q = 0.2. When different initial conditions as (10) and (11) are employed, respectively, the stripe and spot pattern coexist in the space with different structure due to different initial conditions. Of course, we have to mention the fact that due to Hopf bifurcation and Turing instability, the resulting pattern could exhibit oscillations strongly.
In Fig. 10, the control parameter q is set as q = 0.2 with the initial condition of (11). When p takes the value of 0.14, we observe a mixture of cold spots and stripes (cf. Fig. 10a). As p decreases to 0.13, the stripe patterns gradually disappear with the evolution of time t. As a result, cold spots patterns will prevail in the whole domain (cf. Fig. 10b). Finally, we choose the same point in Turing-Hopf region to simulate, that is, choosing q = 0.2 and p = 0.14. While under different initial conditions (11) and uniformly distributed random perturbation separately, the distributions of cold spot and stripe patterns are somewhat different though the snapshots looks alike (cf. Fig. 11a, b).

Conclusions and discussions
To investigate the effect of cross-diffusion on the spatiotemporal patterns of Schnakenberg model, in this paper, we deduce the amplitude equations and present Turing pattern selection using mathematical analysis and numerical simulations. More specifically, it is found that even if d 11 > d 22 , Turing instability can still occur with cross-diffusion. This phenomenon seems negative for systems only with self-diffusion. From the theoretical analysis and the pattern formations through numerical simulations, we can find out that the numerical results are almost in accord with the amplitude equations. But in the circumstance of the stripe pattern formation, the amplitude equations cease to be effective. The main reasons for the invalidation are that μ is far from μ 3 which determines the emergence and the stability of stripe patterns or that the Turing bifurcating parameter q is far away from the critical value q T . We demonstrate that cross-diffusion can lead to spatial and spatiotemporal pattern formation in both Turing instability region and Turing-Hopf region even when the cross-diffusion coefficients are relatively small. In Turing instability region, we show the sequence of pattern transition "cold hexagons → mixture of cold hexagons and stripes → stripes (labyrinthine) → mixture of hot spots and stripes → hot spots" when the Turing bifurcating parameter q ranges between q H and q T . When and what kind of Turing patterns will emerge have been predicted by the amplitude equations. In Turing-Hopf instability region, we also observe the cold spots pattern as well as the mixture of cold spots and stripes patterns. It is worth mentioning that either the uniformly distributed random perturbation or the symmetric initial conditions (10) and (11), for same parameter values, the intrinsic structure of the final patterns stay the same with a little differences in shape. The results obtained in this paper may enrich the pattern dynamics in chemical reaction-diffusion models and help us to better understand the mechanism of crossdiffusion and interactions in a real environment.