The travelling wave of Gray-Scott systems – existence, multiplicity and stability

ABSTRACT This article studies existence and stability of travelling wave of unstirred Gray-Scott system in biological pattern formation which models an isothermal chemical reaction involving two chemical species, a reactant A and an auto-catalyst B, and a linear decay , where C is an inert product. Our result shows a new and very distinctive feature of Gray-Scott type of models in generating rich and structurally different travelling pulses than related models in the literature. In particular, the existence of multiple travelling waves which have distinctive number of local maxima is proved. Furthermore, the stability of travelling wave of the reaction–diffusion system of isothermal diffusion system , is also studied.


Introduction
In this paper, we consider It models chemical reaction of the form A + mB → (m + 1)B with rate uv m and B → C with C being an inert chemical species. Here, D, a positive constant, is the ratio of diffusion coefficients of chemical species B to that of A, m ≥ 1 is a positive constant not necessarily an integer, and kv l describes the rate of B → C, with k and l ≥ 1 both positive constants. We assume throughout that 1 ≤ l < m. Many models in mathematical biology take the form of system (I), see [10,14,17]. In particular, m = 2 and l = 1 is the famous Gray-Scott model in biological pattern formation, one of the popular models proposed for replicating experiment results in early 1990s, see [13,14]. The most exciting feature of the diffusive Gray-Scott system with feeding is self-replicating travelling pulse (travelling wave). It has been extensively studied, see [7,8,11] by formal analysis and numerical computation, but the phenomenon and underlying mechanism is not completely understood. In particular, rigorous analysis of Gray-Scott is badly needed.
Our main concern is the existence and stability of travelling wave to (I) and a related system (II) below. A travelling wave solution for (I) links one equilibrium point to another. Since any equilibrium point of (I) is of the form (a, 0), a ∈ R, the travelling wave problem takes the from u + cu = uv m , u > 0 in R, It is easy to show that the travelling wave solution to (I) exists only when m ≥ l. This is because with u 0 > 0 we need kv m ≥ uv l for all x close to −∞ to make v positive for all x close to −∞. In case of l > m, it is proved in [9] that (I) has a travelling wave solution when u 0 = 0. The dynamics of that model is very different from our case.
The travelling wave problem (1) with m > l, is very different from any of the main types of scalar equation as well as other related models such as the case of l = m, see [6] for more details or the system (II) below of isothermal diffusion system without decay. In particular, for l = m or system (II), the travelling wave problem is of the mono-stable case of scalar equation, but our case is not.
For simplicity, we shall only treat the case of l = 1, m > 1 in details in this work. For the case of m > l > 1, we refer the reader to [20]. By a simple scaling, we can scale out k, and hence, we assume k = 1 here and in Section 2.
The following result is proved in [6] which shows a complete different and unique feature of Gray-Scott type of models.
Theorem A: Let D > 0, m > 1, k = l = 1, and u 0 > 0 be given constants. There exists a positive constant c such that (1) admits a solution. In addition, the set of speeds for existence lies in a bounded interval for a given value of u 0 > 0. Furthermore, the speed c must satisfy That is, our travelling wave problem is not mono-stable type, nor the bistable type as our main result below shows.
Perhaps, the most exciting and surprising result is the one which shows when u 0 1, there exists a large number of travelling wave solutions, each with fixed number of local maxima for w = uv m−1 and with different speed. For this, we re-cast (1), after simple scaling, as where d = D −1 . We make the following change of scale and variables: Define s + = max{s, 0}. Our main result on travelling wave of (I) is the following theorem. (1) There exist positive constants M 1 , M 2 , and M 3 that depend only on m and D such that for each ε > 0, Equation uniformly in i = 1, . . . , L and locally uniformly in z ∈ R, where W is the unique solution of Remark 1.1: It is clear that since u is strictly increasing, the L-hump solution must have at least L local maxima for v, thus a travelling wave solution with a large number of oscillations. As a matter of fact, it should be clear from our proof that as ε → 0, v has exactly L-hump.
In addition, the solution W of Equation (5) is a closed orbit, our result can be explained as that v is a small perturbation of W when ε is small. But, since our system does not admit any close orbit, v has to decay to zero exponentially with a rate uniquely determined by the corresponding travelling wave speed.
Another reaction-diffusion system, we would like to pursue models the following simple isothermal chemical reaction: A + mB → (m + 1)B with rate kab m and m = 1, 2 between two chemical species A and B, where k > 0 is the rate constant. It appears in many chemical wave models of excitable media ranging from the idealized Brusselator to real-world clock reactions such as Belousov-Zhabotinsky reaction, the Briggs-Rauscher reaction, the Bray-Liebhafsky reaction and the iodine clock reaction. In those setting, its importance was recognized pretty early, [10,18].
In this work, we study the travelling wave problem of autocatalytic chemical reaction A + mB → (m + 1)B, which, after simple non-dimensionalization, results in the reaction-diffusion system, For a travelling wave solution to (II), u( where c > 0 is a constant. Assuming It is easy to prove that for a travelling wave solution, lim z→∞ (u, v) = (a, 0). By a simple scaling, we only need to consider a = 1, the travelling wave problem of (II) is the following: It turns out the travelling wave problem of (II) is of the mono-stable type for scalar equation with a minimum positive speed. The main effort is then to prove sharp bound on minimum speed, because the minimum speed travelling wave is the most stable one. The following two results are proved in [4,5]. .
where K(m) is a constant which depends on m only and is an increasing function of m. In particular, K(1) = 1/4, K(2) = 2. For existence, we have the following results.
(i) If m 2, a unique (up to translation) travelling wave solution exists for (7) for each

up to translation) travelling wave solution exists for Equation
Theorem C: Suppose D 1 and m 1. There exists a positive constant c min such that Equation (7) has a travelling wave if and only if c c min . In addition, c min is bounded by

where K(m) is the same constant as in above theorem.
It is clear from the above two results that c min is of order 1. Our main focus on (II) is to study stability of travelling waves using combination of rigorous analysis and numerical computation. Our results on (II) are stated and proved in Section 3.
The organization of the paper is as follows. In Section 2, we study the system (3) and prove the existence of multiple travelling waves to system (I), and in Section 3 we analyse the stability of travelling waves of system (II).

Existence of multiple travelling waves
In this section we study (3) and prove Theorems 1.1. Due to limited space, we present only key steps of the proof and a more detailed presentation will appear later in [15].

Preliminary
For each constant c 0, we consider the initial value problem, for (u, v) = (u(x), v(x)), where λ is the positive root of λ 2 + cελ = 1 and v + := max{v, 0}. The proof is standard and we omit the details. We refer the interested reader to [6]. In the sequel, (u, v) refers to the solution of Equation (8) with (ε, c) dependence suppressed. We use I = (−∞, X) to represent the interval on which v > 0. It is easy to rule out the possibility of u + v tending to ∞ at a finite x. Because if it does, then we must have v tending to ∞ at a finite x. But, this is impossible because of the negative sign of the nonlinear term [ For estimates of v, first we investigate all critical points and possible oscillatory nature of v.
As a consequence, the set {z ∈ R : v (z) = 0, v (z) = 0} can be arranged from small to large by {z i } n i=1 , where either n = ∞ or n is a positive integer. For the latter case, set z 0 = −∞ and z n+1 = ∞. Then for each integer i satisfying Proof: The proof is straight forward.
Next, we establish an explicit upper bound of v.
In particular, taking z = −∞ we have v(x) < M for all x ∈ R, where M is as in Equation (4).

An Upper Bound of c
We now show that there is no travelling wave of fast speed.

Consequently, Equation (3) admits no solution when c max{
Proof: We divide the proof into three steps.

A Differential Inequality.
Let c > 0 be a constant and (u, v) be the unique solution of Equation (8). Set In view of Equation (13), we see that K is finite, so E is well-defined. Using Equations (11) and (12), we derive that Assume that Then we obtain (14). First of all, since 0 < ρ ρ whereρ is the positive root of ε(Dρ 2 + cρ) = M m , we have

A Necessary Condition for Equation
Thus, the first and second inequalities in Equation (14) hold if we have Next, we derive from Equation (13) that Thus, the third inequality in Equation (14) holds if In summary, the assumption (14) holds if c max{ ,

Asymptotic Behaviour as x → ∞.
Now assume that c max{ √ M 1 /ε, M 2 }. Then E + cεE < 0 in R. Consequently, E < 0 in R. Since X < ∞ would imply E(X) > 0, we must have X = ∞, i.e. w > 0 in R. In addition, from E < 0 in R and m > 0, we derive that both w and w are bounded.

Existence of multiple travelling waves for small ε
In the sequel, we always assume that ε and c are parameters satisfying We denote by (u, v) the solution of Equation In addition, setting we have ς = O( √ ε), and Furthermore, exactly one of the following holds:

â). In this caseâ is a local positive minimum of w;
(ii) w(â) = 0 > w (â). In this case,â = X < ∞ and w < 0 in [X, ∞); For detailed proof, we refer the reader to [15]. We define an energy functional by Note that E = w 2 − G(w) + O(ε)w 2 and when 0 w 1, Direct differentiation together with the differential equation for w gives Let a 1 = −∞ be the 'first local minimum' of w and b 1 be the first local maximum. We set r 1 = w(a 1 ) = w(−∞) = 0 and R 1 = w(b 1 ). Then by Lemma 2.6, w > 0 in (a 1 , b 1 ) and R 1 = M + O(ς ).
Let i = O(1)/ √ ε be a positive integer and assume that w has at least i local minima attained, from small to large, at a 1 , a 2 , . . . , a i , satisfying r j := w(a j ) ∈ [0, 1/2] for j = 1, . . . , i. Then by Lemma 2.6, there exist i local maxima, attained, from small to large, at a i+1 ) is the maximal interval on which w < 0 < w. Set r i+1 = w(a i+1 ). By Lemma 2.6, we have G(r i+1 ) = G(r i ) + O(ς i ). We call γ i := {(w(x), w (x)) : x ∈ (a i , a i+1 )} the ith loop of the trajectory on the w − w phase plane. We observe the following: E(a i+1 ) = 0, then w(a i+1 ) = 0 and w (a i+1 ) = 0 so we have a solution of Equation (8) Hence, w(a i+1 ) is a local minimum of w and the trajectory has at least i+1 loops on the phase plane.
In the sequel, we evaluate E(a i+1 ) in terms of E(a i ). For each positive integer n not too large, we shall find an appropriate c = c n > 0 such that E(a i+1 ) < 0 for 1 i < n and E(a n+1 ) = 0, i.e. the solution of Equation (8) is a solution of Equation (3) with exactly n loops; as a consequence, we obtain an n-hump travelling wave.
Integrating Equation (18) over (a i , a i+1 ), we find that a i w 2 (y) dy, Through very tedious computation, we can prove the following: ) and (b i , a i+1 ) be the maximum interval on which w < 0 < w. Then

Lemma 2.7: Suppose a i is the ith local minimum of w, w(a i ) ∈ [0, 1/2], and b i is the ith local maximum. Let r i = w(a i
We now evaluate ρ(b 1 ) and the minimal speed wave. Note that a 1 = −∞, r 1 = 0, and c).
Thus, we have the following: (4). Then

Lemma 2.8: Let σ and γ be as in Equation
Now we are ready to prove the following:

The stability of travelling waves
In this section, we study the stability of travelling wave solutions to (II) in R.

Analytical results
In spite of the apparent importance and close relation to the classical Fisher-KPP scalar equation, there are very few analytical results on (II) with m > 1. For some of the most recent development, see [5,12]. For convenience, we do a change of variables, and after dropping ' ∼ ', we have The very important issue for (II), in light of the experiment [18], is how fast the spreading of local disturbance of v under the laboratory initial conditions of u 0 (x) ≡ 1 and v 0 (x) has a compact support. Our main analytical result is the following theorem.

The computational approach
In this part, we present some computational results on (II) with two special cases m = 1 and m = 2, and the diffusion coefficient D either in (0, 1) or in (1, ∞). The purpose is two folds. On the one hand, computation can verify and confirm analytical results, in particular, whether the spreading of local disturbance of v is of order O( √ D) when D > 1. On the other hand, it can help us to gain insight into the complex interaction of diffusion and nonlinear reaction terms and how their interaction determines the behaviour of solutions. This is very important in our study of the stability of travelling waves and the limiting profiles of solution as t → ∞.
In all examples of computation, the initial conditions are u(x, 0) = 1 and v(x, 0) has a compact support. We take the spatial domain to be a large interval centred at zero and use periodic boundary conditions. Figure 1 is the result of computation of m = 1, D = 2 with initial condition of v to be v(x, 0) = 1 in [−1, 1], and zero otherwise. The spatial domain is [−40, 40]. The reaction starts from the central region and spread out with the speed c, which is approximately 2 √ D, the minimum speed, before v becomes very flat, approaching 1. This is in agreement with the theoretical result of [5]. Figure 2 is the result of computation of m = 2 with the other conditions same as the above case. The reaction again starts from the central region and spread out with the estimated speed of 2.5 √ D/3, before v becomes very flat, approaching 1 as time t 1. Figure 3 is the result of computation of m = 2 and D = 4 with other conditions same as the above case. The reaction again starts from the central region and spread out with the estimated speed of 2.5 √ D/3, before v becomes very flat, approaching 1 as time t 1.        converges to a fixed bell-shaped profile as time t 1. In addition, u becomes two-hump from the initial one-hump and keep the same profile with diminished height as t increases before eventually tending to zero. Figure 5 is the result of computation of m = 2 and D = 3 with other conditions same as the above case. The solutions demonstrate the same kind of qualitative behaviour as the above case except the speed range is in 7 √ D/12 < c < 5 √ D/6. We also did some computation on (I) with the results shown in the figures below. In Figures 6-8, we present some computational results on (I) with m = 2, l = 1, D = 4 and the same initial conditions as for various cases of (II). They show that when k = 1 and k = 0.2, the decay is very strong and v tends to zero very fast before any pattern to form effectively. But, for k = 0.05, v undergoes some very interesting evolution before decay to zero eventually.