Stable spike clusters for the one-dimensional Gierer–Meinhardt system

We consider the Gierer–Meinhardt system with precursor inhomogeneity and two small diffusivities in an interval $$\begin{equation*}
\left\{
\begin{array}{ll}
A_t=\epsilon^2 A''- \mu(x) A+\frac{A^2}{H}, &x\in(-1, 1),\,t>0,\\[3mm]
\tau H_t=D H'' -H+ A^2, & x\in (-1, 1),\,t>0,\\[3mm]
A' (-1)= A' (1)= H' (-1) = H' (1) =0,
\end{array}
\right.
\end{equation*}$$ $$\begin{equation*}\mbox{where } \quad 0<\epsilon \ll\sqrt{D}\ll 1, \quad
\end{equation*}$$ $$\begin{equation*}
\tau\geq 0
\mbox{ and
$\tau$ is independent of $\epsilon$.
}
\end{equation*}$$ A spike cluster is the combination of several spikes which all approach the same point in the singular limit. We rigorously prove the existence of a steady-state spike cluster consisting of N spikes near a non-degenerate local minimum point t 0 of the smooth positive inhomogeneity μ(x), i.e. we assume that μ′(t 0) = 0, μ″(t 0) > 0 and we have μ(t 0) > 0. Here, N is an arbitrary positive integer. Further, we show that this solution is linearly stable. We explicitly compute all eigenvalues, both large (of order O(1)) and small (of order o(1)). The main features of studying the Gierer–Meinhardt system in this setting are as follows: (i) it is biologically relevant since it models a hierarchical process (pattern formation of small-scale structures induced by a pre-existing large-scale inhomogeneity); (ii) it contains three different spatial scales two of which are small: the O(1) scale of the precursor inhomogeneity μ(x), the $O(\sqrt{D})$ scale of the inhibitor diffusivity and the O(ε) scale of the activator diffusivity; (iii) the expressions can be made explicit and often have a particularly simple form.


Introduction
In his pioneering work [31] in 1952, Turing studied how pattern formation could start from an unpatterned state. He explained the onset of pattern formation by the presence τ > 0 and τ is independent of ε.
In the standard Gierer-Meinhardt system without precursor, it is assumed that μ(x) ≡ 1.
Precursor gradients in reaction-diffusion systems have been investigated in earlier work. The original Gierer-Meinhardt system [11,16,18] has been introduced with precursor gradients. These precursors were proposed to model the localization of the head structure in the coelenterate Hydra. Gradients have also been used in the Brusselator model to restrict pattern formation to some fraction of the spatial domain [14]. In that example, the gradient carries the system in and out of the pattern-forming part of the parameter range (across the Turing bifurcation), thus effectively confining the domain where peak formation can occur.
In this paper, we study the Gierer-Meinhardt system with precursor and prove the existence and stability of a cluster, which consists of N spikes approaching the same limiting point.
More precisely, we prove the existence of a steady-state spike cluster consisting of N spikes near a non-degenerate local minimum point t 0 of the positive inhomogeneity μ(x) ∈ C 3 (Ω), i.e. we assume that μ (t 0 ) = 0, μ (t 0 ) > 0 and we have μ(t 0 ) > 0. Further, we show that this solution is linearly stable.
We explicitly compute all eigenvalues, both large (of order O (1)) and small (of order o (1)). The main features of studying the Gierer-Meinhardt system in this setting are as follows: (i) It is biologically relevant since it models a hierarchical process (pattern formation of small-scale structures induced by a pre-existing inhomogeneity).
(ii) It is important to note that this system contains three different spatial scales two of which are small (i.e. o (1) To have a non-trivial spike cluster, we assume throughout the paper that Before stating our main results, let us review some previous results on pattern formation for the Gierer-Meinhardt system (1.1) , in particular concerning spiky patterns.
(1) I. Takagi [30] proved the existence of N-spike steady-state solutions of (1.1) in an interval for homogeneous coefficients (i.e. μ(x) = 1) in the regime ε 1 and D 1, where N is an arbitrary positive integer. For these solutions, the spikes are identical copies of each other and their maxima are located at The proof in [30] is based on symmetry and the implicit function theorem.
(2) In [15] (using matched asymptotic expansions) and [43] (based on rigorous proofs), the following stability result has been shown: For the N-spike steady-state solution derived in item 1 and 0 6 τ < τ 0 (N), where τ 0 (N) > 0 is independent of ε, there are numbers D 1 > D 2 > · · · > D N > · · · (which have been computed explicitly) such that the N-spike steady state is stable for for D < D N and unstable for D > D N .
(3) In [15] (using matched asymptotic expansions) and [33] (based on rigorous analysis), the following existence and stability results have been shown: For a certain parameter range of D, the Gierer-Meinhardt system (1.1) with μ(x) = 1 has asymmetric N-spike steady-state solutions, which consist of exact copies of precisely two different spikes with distinct amplitudes. They can be considered as bifurcating solutions from those in item 1 such that the amplitudes start to differ at the bifurcation point (saddle-node bifuraction). The stability of these asymmetric N-peaked solutions has been studied in [33].
(4) In [45], the existence and stability of N-peaked steady states for the Gierer-Meinhardt system with precursor inhomogeneity has been shown. These spikes have different amplitudes. In particular, the results imply that precursor inhomogeneities can induce instabilities. Single-spike solutions for the Gierer-Meinhardt system with precursor including spike motion have been studied in [32].
(5) In [42], the existence of symmetric and asymmetric multiple spike clusters in an interval has been shown.
Compared to each of the items listed above, the setting and results in our paper have marked differences. We now consider two small parameters, D and ε √ D which results in new types of behaviour. The leading-order asymptotic expression of the large and small eigenvalues depend on the index of the eigenvalue quadratically, whereas in items 1 and 2, this relation is oscillatory (involving trigonometric functions). In our study, the spikes in leading order have equal amplitudes and uniform spacing, although there is precursor inhomogeneity in the system, in contrast to item 3. The amplitudes, positions and eigenvalues in our study can be characterized explicitly and have a simpler form than in item 4. We can also prove the stability of clusters not merely their existence as in item 5. In particular, we show here that the clusters may be stable, whereas in item 5 they are expected to be unstable. In the shadow system case (D = ∞), the existence of single-or N-peaked solutions has been established in [12,13,21,22] and other papers. It is interesting to remark that symmetric and asymmetric patterns can also be obtained for the Gierer-Meinhardt system on the real line, see [7,8]. We refer to [23] for the singular limit eigenvalue problem (SLEP) approch for the existence and stability of multi-layered solutions for reaction-diffusion systems. For two-dimensional domains, the existence and stability of multi-peaked steady states has been proved in [38][39][40] and results similar to items 1 and 2 have been derived. Hopf bifurcation has been established in [6,34,35,40]. The repulsive dynamics of multiple spikes has been studied in [9].
Another study with three different spatial scales, two of which are small, considers a consumer chain model allowing for a novel type of spiky clustered pattern which is stable for certain parameters [46].
The model in our paper shows some similarity to variational models for material microstructure [1,20,48]. In both models, the solutions have two small scales. However, in our case, we have two parameters to control each of them independently, whereas in the microstructure case they are expressions of different orders depending on the same small parameter and so they are related to one another.
Results on the existence and stability of multi-spike steady states have been reviewed and put in a general context in [47].
We plan to continue the investigation of stable clusters in future research. In particular, we are interested in two-dimensional patterns. Whereas in one space dimension the spikes in the cluster are aligned, in two dimensions we expect a rather rich geometric picture of possible spike locations. This paper has the following structure: In Section 2, we state our main results on existence and stability. Then, we show some numerical simulations to illustrate the main results. We also study the dynamics of pattern formation, even outside the regime covered by the results. Next, we present four highlights of the results and their proofs. Finally, we sketch the main steps of the proofs. In Section 3, we introduce some preliminaries. In Sections 4-7, we prove the existence of steady-state spike clusters: In Section 4, we introduce suitable approximate solutions, in Section 5, we compute their error, in Section 6, we use the Liapunov-Schmidt method to reduce the existence of solutions of (1.2) to a finite-dimensional problem, in Section 7, we solve this finite-dimensional reduced problem. In Sections 8, 9 and Appendix B, we prove the stability of these steady-state spike clusters: In Section 8, we study the large eigenvalues of the linearized operator and show that it has diagonal form. We give a complete description of their asymptotic behaviour which is stated in Lemma 15. In Section 9, we characterize the small eigenvalues of the linearized operator and show that they all have negative real part. This includes deriving the eigenvalues of a matrix which is needed to compute the small eigenvalues explicitly. We give a complete description of their asymptotic behaviour in leading order which can be found in Lemma 16. Our approach here is to interpret the main matrix as the finitedifference approximation of a suitable ordinary differential equation, compute the solution of this approximation explicitly and get the eigenvectors by taking its values at uniformly spaced points. In Section 10, we conclude with a discussion of our results with respect to the bridging of length scales and the hierarchy of multi-stage biological processes. In Appendix A, we state a few results on NLEPs which are needed throughout the paper. In Appendix B, we perform the technical analysis needed to derive the small eigenvalues.

Main results on existence and stability
In this section, we state our main results on existence and stability of solutions and present four highlights of our approach and sketch the proofs of the main results.
Let t 0 ∈ (−1, 1) and set compare (4.30). We set Our first result is about the existence of an N-spike cluster solution near a non-degenerate local minimum point of the precursor.
Next, we state our second result which concerns the stability of the N-spike cluster steady states given in Theorem 1.
Theorem 2 (Stability of an N-spike cluster) For ε √ D 1, let (A ε , H ε ) be an N-spike cluster steady state given in Theorem 1. Then, there exists τ 0 > 0 independent of ε and √ D such that the N-spike cluster steady state (A ε , H ε ) is linearly stable for all 0 6 τ < τ 0 .  Figure 1. Clustered spiky steady states of (1.1) for ε 2 = 0.00001, D = 0.001, μ(x) = 6 + 400x 2 and μ(x) = 6 + 200x 2 , respectively. Shown are a six-spike cluster and an eight-spike cluster. In both cases, the activator a is displayed in the left graph and the inhibitor h in the right graph. Now by performing some numerical simulations, we study patterns for the Gierer-Meinhardt system with precursor gradient given in (1.1) systematically in various situations.
Throughout all these simulations, we choose τ = 0.1, ε 2 = 0.00001 and take varying values for D ranging between a few times ε 2 and 1.
Further, we vary the strength of the precursor gradient and observe two distinct types of behaviour: For strong precursor gradient, the spikes assemble as a cluster near the global minimum point of the precursor, for weak precursor gradient, the spikes are distributed over the whole interval.
We will observe a rather rich dynamical behaviour which by far exceeds the immediate vicinity of the spike cluster which will be analysed in detail in this paper.
First, we show the results of computations of spike cluster steady states of (1.1). These have been obtained as long-time limits of (1.1) and are numerically stable.
In Figure 1, we note that the amplitudes of the different spikes of a cluster in the activator component are very close to each other.
On the other hand, the inhibitor values differ substantially at different spike locations and are highest near the centre of the cluster. This stands in contrast to the precursor μ(x) which attains its global minimum at the centre of the cluster. The combination of these two effects leads to almost equal spike amplitudes.
Both the activator and inhibitor peaks have almost equal distance. Now, we show the initial condition which has been used in all simulations of clusters and multiple spikes.
Next, we display the dynamics of getting a cluster from the initial condition given in Figure 2.
In Figure 3, for t = 0, we start with very fine oscillations of the activator a. Then, for t = 0.70, we reach a pattern which is very close to zero except in the centre of the interval where fine oscillations still prevail. Starting from t = 0.80, we see eight spikes which generally increase in amplitude and whose positions are mainly fixed. Looking more closely, we can also see that the amplitudes show some oscillatory behaviour, first overshooting the final amplitude of around 0.57 and then oscillating around the final amplitude and approaching it.
Next, we consider a different regime of a weaker precursor gradient (all other parameters remain unchanged from before). We observe spikes which are distributed over the whole interval.  Now, we display the dynamics of getting a multiple spike pattern from the initial condition given in Figure 2.
In Figure 5, we start with very fine oscillations of the activator a. Immediately, the amplitudes of the oscillations change driven by the precursor gradient, but their period remains mainly unchanged. Then, spikes start to form, first at the boundary, then further and further inside the interval. Finally, the spikes in the centre of the domain form almost simultaneously and there is also some oscillatory behaviour between their amplitudes. It takes longest for the amplitudes of the two spikes in the centre to increase to their steady-state value. The positions of the spikes are mainly unchanged but the increase in amplitude for the two central spikes is coupled with them slowly moving apart and pushing the remaining spikes away from the centre.
In Figure 6, we show how the number of multiple spikes depends on the diffusion constants.
Finally, in Figure 7, we show the behaviour as the two small diffusion constants ε 2 and D come close to each other.

Remark 3
For the stability, we assume that 0 6 τ < τ 0 for some τ 0 > 0. Stability in the case where τ is large has been investigated in [35] for single spikes and those results on Hopf bifurcation are expected to carry over to the case of an N-spike cluster considered here. We remark that stability in the case of large τ for the shadow system has been studied in [6,34]. It turns out that this Hopf bifurcation leads to oscillations of the amplitudes. The Hopf bifurcation at τ = τ 0 , where τ 0 is of order 1, still arises even in the regime 0 ε √ D 1 considered in this paper. For the spike cluster, τ 0 is independent of N which can be shown by the analysis in Sections 8 and 9.

Remark 4
It is an interesting question to consider the maximum number of stable spikes in the regime 0 ε √ D 1 studied in this paper. We expect that there are stable multispike solutions if N < c √ D (in leading order of D). In the regime N ∼ c √ D , the spikes will be distributed over the whole interval and at N ∼ c √ D , we expect an overcrowding instability. This threshold would be an extension of the corresponding result for D = 1 (see [43]). The cluster solutions studied in this paper are only possible if N < c √ D log 1 D (in leading order of D) due to the distance between spikes. We have presented some numerical simulations of

Remark 5 If ε
√ D ∼ 1, the spike solutions will change into other types of patterns, e.g. spatial oscillations, which could again be stable. It is also possible that the patterns will vanish. We have presented some numerical examples to illustrate this behaviour (see Figure 7).

Remark 6
Previous studies of the precursor case can be found amongst others in [2,27,28]. We also refer to results for the dynamics of pulses in heterogeneous media [24,49]. This clustered spike pattern and multiple spikes distributed over the whole interval are more regular than multiple spike patterns observed for the Gierer-Meinhardt system with precursor  and order 1 inhibitor diffusivity studied in [45]. In that case, multiple spike patterns have irregular distances and amplitudes since the precursor interacts with the geometry of the domain (represented by the Green's function) globally. On the other hand, for the regime covered in this paper, the precursor acts globally but the Green's function acts only locally between neighbouring spikes.
The proofs of both Theorems 1 and 2 will follow the approach in [47], where we reviewed and discussed many results on the existence and stability of multi-spike steady states. Before providing a sketch of the proofs of Theorems 1 and 2, we first state four highlights of the results and proofs in an informal manner. For each of these, we also indicate the novelty in comparison to previous work.
Highlight 1: For the proof in Theorem 1, we use Liapunov-Schmidt reduction to derive a reduced problem which will determine the positions of the spikes. This reduced problem in leading order is given by where c 1 , c 2 > 0 are constants which are independent of the small parameters and which implies (compare (7.6)). The distance between neighbouring spikes in the cluster is small (converging to zero) and in leading order, it is the same between any pair of neighbours. Further note that the distance is determined by the large diffusion constant D. In previous work for D = O(1), the spikes either have distance of order 1 or, in the case of clusters, they have small distance which is determined by the small diffusion constant ε. Thus, the spacing between neighbouring spikes follows a new asymptotic rule not encountered before.
The deeper reason for this behaviour is that the spike cluster we consider in this paper is formed by balancing the interactions of the inhibitor between neighbouring spikes and the inhomogeneity. On the other hand, in [42] and other previous work on spike clusters, they are established by balancing the interactions of the activator between neighbouring spikes and the inhomogeneity.
Highlight 2: The large eigenvalues with λ ε → λ 0 = 0 and their corresponding eigenfunctions where φ ε,i (y) is the restriction of the re-scaled eigenfunction of the activator A ε near t i , in (see (8.6)). This NLEP has diagonal form. Thus, with respect to large eigenvalues each spike only interacts with itself and not with the other spikes. It follows that the spike cluster is stable with the respect to large eigenvalues. In previous work for the case D = O(1), the stability problem of large eigenvalues for multiple spikes leads to a vectorial NLEP. It has to be studied by considering the spectrum of the matrix involved. Depending on the parameters, the multiple spikes can be stable or unstable. In the case of clusters for D = O(1), the stability has not been considered rigorously but we expect that the solution is unstable. Highlight 3: The small eigenvalues λ ε → 0 in leading order are given by the eigenvalues of the matrix where c 3 > 0 is independent of the small parameters, t 0 = (t 0 , . . . , t 0 ) and with δ N,0 = δ 1,N+1 = 0 (compare 9.13). The tridiagonal matrix M(t 0 ) derived here indicates that with respect to small eigenvalues each spike only interacts with its nearest neighbour and not with the other spikes. This is different from the case D = O(1): For multiple spikes, the matrix typically has strictly positive entries only (although it could have zero eigenvalues in the presence of symmetries), for clusters, a similar tridiagonal matrix which depends on ε has been studied to show the existence of spike clusters, see [42].

Highlight 4:
We determine all the eigenvalues of the matrix M(t 0 ) (see Highlight 3) explicitly by a method based on exactly finding a finite-difference approximation to a suitable ordinary differential equation.
These eigenvalues are given by Further, there is an eigenvalue of smaller size given by Lemma 16). This implies that the spike cluster is stable with respect to small eigenvalues. This result seems to be new in the literature. Therefore, we provided this proof.
Finally, we give a sketch of the proofs of Theorem 1 (existence) and Theorem 2 (stability).
We begin by stating the main steps in the proof of Theorem 1: (1) The existence problem (4.20) is reduced to a non-local one-dimensional problem is an integral operator solving (4.21).
(2) The ansatz for a spike cluster isÂ (3) The amplitude identity is crucial to show that S ε [Â] is small in an appropriate norm. Therefore,Â is intuitively almost a solution.
(4) Using the estimate on S ε [Â], one can perform Liapunov-Schmidt reduction resulting in a small such that ⊥ C amounts to equating N L 2 -inner products to zero. This leads to a system of equations to leading order given by (2.11) in the spike points (t 1 , . . . , t N ), which can be shown to have a solution.
The main steps in the proof of Theorem 2 can be stated as follows: The eigenvalue problem is derived by linearizing the Gierer-Meinhardt system (1.1) around the clustered steady state A ε derived in Theorem 1. It is stated in (8.1) and for τ = 0, we havẽ where λ ε is some complex number and φ ε ∈ H 2 (Ω) satisfying Neumann boundary conditions. Further, for φ ∈ L 2 (Ω), is an integral operator solving (8.3). Then, we consider the eigenvalues in three cases separately as follows: (1) We first study large eigenvalues λ ε = O(1).
(i) Using (8.1) and the decay of the spikes, it is shown that in leading order an eigenfunction satisfies This means that the eigenfunction can be decomposed into parts which are located near each of the spikes.
This means that we get an NLEP which has diagonal form.
(iii) By a result in [43], it follows that the spike cluster is stable with respect to large eigenvalues.
(i) The ansatz for an eigenfunction is (ii) Then, we show that φ ⊥ = o (1). Note that εw ε,j is of exact order 1 and so it dominates the eigenfunction.
(iii) Taking the spatial derivative of the steady-state problem (4.20), we get an identity forw ε,j (up to small error terms) which we subtract from the eigenvalue problem.
(iv) We expand the terms in (8.1) around the spikes using the expansion of the Green's function G D and the inhomogeneity μ(x), both around the spike points, collecting the leading terms and giving rigorous estimates for the remainder.
Using the results in (iv), we derive the matrix M(t 0 ) stated in (9.13) which determines the stability properties caused by small eigenvalues.
(vi) We determine the eigenvalues and eigenvectors of the matrix M(t 0 ) explicitly by considering it as the finite-difference approximation of a suitable ordinary differential equation, compute the solution of this approximation explicitly and get the eigenvectors by taking its values at uniformly spaced points. See Lemma 13 and its proof.
It follows that the eigenvalue problem is stable with respect to small eigenvalues.
(3) Finally, we show that there are no eigenvalues λ ε with |λ ε | → ∞, if we assume that Re(λ ε ) > −c for some fixed c > 0. (If the latest condition is not satisfied the eigenfunction will be uniformly bounded in time and so does not cause any instability.) (i) We multiply (8.1) on both sides by φ ε . This leads to a quadratic form in φ ε .

Preliminaries: scaling property, Green's function and eigenvalue problems
In this section, we will provide some preliminaries which will be needed later for the existence and stability proofs. Let w be the ground state solution given in (2.1). By a simple scaling argument, the function is the unique solution of the problem z) (and the left-hand and right-hand limits are considered for x = z). We calculate to be the singular part of G D (x, z). Let the regular part H D of G D be defined by where d 0 = 1 − |t 0 | and η 0 > 0 is an arbitrary but fixed constant. Forξ 0 , we estimatê Let us denote ∂ ∂ti as ∇ ti . When i = j, we can define ∇ ti G(t i , t j ) in the classical way because the function is smooth.
Similarly, we define (3.9) For convenience and clarity, we introduce a re-scaled version of the Green's function which has a finite limit as D → 0. Thus, we set Throughout the paper, let C, c denote generic constants which may change from line to line.

Existence proof I: approximate solutions
Let t 0 ∈ (−1, 1) be a non-degenerate local minimum point of the precursor inhomogeneity, i.e. we assume that (2.7) is satisfied. In this section, we construct an approximation to a spike cluster solution to (1.2) which concentrates at t 0 .
The approximate cluster consists of the spikes μ i w √ μ i x−ti ε which are centred at the points t i and have the scaling factors and η > 0 is a constant which is small enough and will be chosen in Section 7 (see equation (7.9)). The reason for assuming (4.2) and (4.3) will become clear in Section 7 when we solve the reduced problem. We further denote and set To simplify our notation, for t ∈ Ω η and k = 1, . . . , N, we set where χ is a smooth cut-off function which satisfies the conditions and ε δ ε 20 Further, we compute, using (3.5), where d 0 = min(1 − t 0 , t 0 + 1), η 0 > 0 is an arbitrary but fixed constant (compare (3.7)). We haveĜ Generally, we havê For the derivatives, we estimate and analogous results hold for the mixed derivatives. By rescalingÂ = ξ ε A,Ĥ = ξ ε H with ξ ε defined in (2.6), it follows that (1.2) is equivalent to the following system for the re-scaled functionsÂ,Ĥ: (4.20) From now on, we shall work with (4.20) and drop the hats. Next, we rewrite (4.20) as a single equation with a non-local term.
For a function A ∈ H 2 (−1, 1), we define T [A] to be the solution of For t ∈ Ω η , we define an approximate solution to (4.22) by the ansatz where t ∈ Ω η ,w k has been defined in (4.7) andξ k satisfies the amplitude identitŷ Intuitively, this ansatz is close to a solution of (4.22) since by the choice ofw k , the first equation is approximately satisfied and due to (4.24) the second equation holds approximately.
Next, we are now going to determine the amplitudesξ k to leading order so that (4.24) will be satisfied. From (4.21), we have We have O(ε 10 )).
In the section, we will compute its error.

Existence proof II -error of approximate solution
In this section, we compute the error terms caused by the approximate solutions in Section 4. We begin by considering the spatial dependence of the inhibitor near the spikes which is given by . . . , x N ) ∈ Ω η and t ∈ Ω η , where the non-local operator T [A] has been defined in (4.21) and the approximate solution has been introduced in (4.23).
To simplify our notation, we let where ξ ε has been introduced in (2.6) and J 1 is defined by the latest equality. For J 1 , we have by (4.7). We further compute using (4.14). Let x = t k + εz. For k = s, we have where w μs has been defined in (3.1) and is an even function, using (2.6) and (4.10). For k = s, we have using (4.10) and (2.6). Combining (5.5) and (5.7), we have

Remark 7
(i) The second line in (5.8) is an even function in the inner variable y which will drop out in many subsequent computations due to symmetry.
(ii) The third line in (5.8) is an odd function in the inner variable y. For t ∈ Ω η , we have Thus, the third line in (5.8) is of exact order O ε √ D log 1 D y .
Next, we compute and estimate the error terms of the Gierer-Meinhardt system (4.20) for the approximate solution w ε,t . We recall that a steady state for (4.20) is given by and T [A] is defined by (4.21), combined with Neumann boundary conditions A (−1) = A (1) = 0. We now compute the error term Using the amplitude identityξ k = T [A](t k ) and the equation (3.2) of the spike profile for a = μ s , we have Using μ (t 0 ) = 0, we get Now, we readily have the estimate (5.11)

Remark 8
The estimates derived in this section will be needed to conclude the existence proof using Liapunov-Schmidt reduction in Section 6. In particular, they will imply an explicit formula for the positions of the spikes in Section 7.

Existence proof III -Liapunov-Schmidt reduction
In this section, we study the linear operator defined bỹ where A = w ε,t and T [A] has been defined in (8.3).
We will prove results on its invertibility after suitable projections. This will have important implications on the existence of solutions of the non-linear problem including bounds in suitable norms. The proof uses the method of Liapunov-Schmidt reduction which was also considered in [10,12,13,25,26,37] and other works.
We define the approximate kernel and co-kernel of the operatorL ε,t , respectively, as follows: Recall that the vectorial linear operator L has been introduced in (A 2) as follows: By Lemma 21, we know that with X 0 = span{ dw dy } is invertible and possesses a bounded inverse. We also introduce the orthogonal projection π ⊥ ε,t : L 2 (Ω) → C ⊥ ε,t and study the operator L ε,t := π ⊥ ε,t •L ε,t , where orthogonality has been defined in L 2 (Ω) sense. We will show that L ε,t : K ⊥ ε,t → C ⊥ ε,t is invertible with a bounded inverse provided max( ε √ D , D) is small enough. In proving this, we will use the fact that this system is the limit of the operator L ε,t as max( ε √ D , D) → 0. This statement is contained in the following proposition.

Proposition 9
There exist positive constantsδ, λ such that for max ε √ D , D ∈ (0,δ) and all t ∈ Ω η , we have Further, the map We define φ ε,i , i = 1, 2, . . . , N and φ ε,N+1 as follows: At first (after rescaling), the functions φ ε,i are only defined on Ω ε . However, by a standard result, they can be extended to R such that their norm in H 2 (R) is bounded by a constant independent of ε, D and t for max( ε √ D , D) small enough. In the following, we will study this extension. For simplicity of notation, we keep the same notation for the extension. Since for i = 1, 2, . . . , N, each sequence {φ k i } := {φ ε k ,i } (k = 1, 2, . . .) is bounded in H 2 loc (R) it has a weak limit in H 2 loc (R), and therefore also a strong limit in L 2 loc (R) and L ∞ loc (R). Call these limits φ i . Then, passing to the limit in the equation (6.5) in each of the sets Ω j = {x ∈ Ω : x = t j + εy, |y| 6 δε 2ε } (we refer to Appendix A of [40] for further details), Therefore, we conclude φ N+1 = 0 and φ k N+1 H 2 (R) → 0 as k → ∞. This contradicts φ k H 2 (Ωε k ) = 1. To complete the proof of Proposition 9, we just need to show that the conjugate operator to L ε,t (denoted by L * ε,t ) is injective from K ⊥ ε,t to C ⊥ ε,t . The proof for L * ε,t follows along the same lines as for L ε,t and is omitted. Now, we are in the position to solve the equation , we can rewrite this as where where λ > 0 is independent of δ > 0, ε > 0, D > 0 and c(δ) → 0 as δ → 0. Similarly, we show that where c(δ) → 0 as δ → 0. If we choose then for max √ D , D small enough, the operator M ,t is a contraction on B ε,δ . The existence of a fixed point φ ,t now follows from the standard contraction mapping principle and φ ,t is a solution of (6.9).
We have thus proved Lemma 10 There exists δ > 0 such that for every pair of , t with 0 < < δ and t ∈ Ω η , there exists a unique φ ,t ∈ K ⊥ ,t satisfying S [w ε,t + φ ε,t ] ∈ C ε,t . Furthermore, we have the estimate Following the same decomposition into leading even and odd terms as discussed in Remark 7 (see also (5.10)) and applying the linear operator L ε,t to both of them, we get where φ ε,t,1 is an even function in the inner variable y which can be estimated as and φ ε,t,2 is an odd function in the inner variable y which can be estimated as Note that the even term is bigger than the odd term but it will drop in many subsequent calculations.

Existence proof IV: reduced problem
In this section, we solve the reduced problem. This completes the proof of our main existence result given by Theorem 1. By Lemma 10, for every t ∈ Ω η , there exists a unique solution φ ,t ∈ K ⊥ ,t such that To this end, let Then, the map W (t) is continuous in t ∈ Ω η and it remains to find a zero of the vector field W ε (t).
We compute We first compute the main term given by Let x = t s + εy. By (5.10), we can decompose S ε [w ε,t ] into odd and even functions. In leading order, only the odd components of S ε [w ε,t ] matter and we have In summary, we have Next, we estimate Integration by parts and and using the derivative of (3.2) gives This implies which follows from (6.11) and (6.12). Now, for given small ε > 0, we have to determine t ε ∈ Ω η such that W ε,s (t ε ) = 0 for s = 1, . . . , N.
We first consider the limiting case which only takes into account the leading terms and set By (7.6) and (7.7), we have t * ∈ Ω η if D is small enough. We need to find t ε ∈ Ω η such that W ε (t ε ) = 0. Setting e = (1, 1 . . . , 1) T , we have For t ∈ Ω ε , we expand Noting that for t * + τ ∈ Ω η , we have |v| 6 η √ D and α 6 η √ D log 1 D , we get This implies Setting τ = t − t * , we have to determine τ such that and so τ must be a fixed point of the mapping Using projections, we have We now determine when the mapping M ε,D maps from B 1 into B 1 for max( ε √ D , D) small enough. We need to have and we assume By Brouwer's fixed point theorem, the mapping M ε,D possesses a fixed point τ ε ∈ B 1 . Then, t ε = t * + τ ε ∈ Ω η is the desired solution which satisfies W ε (t ε ) = 0.
Thus, we have proved the following proposition.

Proposition 11
For max √ D , D small enough, there exist points t ∈ Ω η with t → t 0 such that W (t ε ) = 0.
Finally, we complete the proof of Theorem 1.

Stability proof I: large eigenvalues
In this section, we study the large eigenvalues which satisfy λ ε → λ 0 = 0 as max( ε Then, we need to analyse the eigenvalue problem where λ ε is some complex number, A ε = w ε,t ε +φ ε,t ε with t ε ∈ Ω η determined in the previous section, and for φ ∈ L 2 (Ω), the function T [A]φ is defined as the unique solution of First, we consider the special case τ = 0. Because we study the large eigenvalues, there exists some small c > 0 such that |λ ε | > c > 0 for max( √ D , D) small enough. We are looking for a condition under which Re(λ ε ) 6 c < 0 for all eigenvalues λ ε of (8.1), (8.2) if max( ε √ D , D) is small enough, where c is independent of ε and D. If Re(λ ε ) 6 −c, then λ ε is a stable large eigenvalue. Therefore, for the rest of this section, we assume that Re(λ ε ) > −c and study the stability properties of such eigenvalues.
Proof (1) of Theorem 12 follows by asymptotic analysis similar to Section 6.
To prove (2) of Theorem 12, we follow a compactness argument of Dancer [6]. The main idea of his approach is as follows: Let λ 0 = 0 be an eigenvalue of problem (8.6) with Re(λ 0 ) > 0.
Then, we can rewrite (8.1) as follows: where R ε (λ ε ) is the inverse of −Δ + (μ(x) + λ ε ) in H 2 (R) (which exists if Re(λ ε ) > − min x∈R μ(x) or Im(λ ε ) = 0) and the non-local operators have been defined in (4.21) and (8.3), respectively. The main property is that R ε (λ ε ) is a compact operator if max( ε √ D , D) is small enough. The rest of the argument follows in the same way as in [6].
We now study the stability of (8.1), (8.2) for large eigenvalues explicitly and prove Theorem 2.
When studying the case τ > 0, we have to deal with NLEPs as in (A 1), for which the coefficient γ of the non-local term is a function of τα. Let γ = γ(τα) be a complex function of τα. Let us suppose that γ(0) ∈ R, |γ(τα)| 6 C for Re(α) = α R > 0, τ > 0, (8.8) where C is a generic constant which is independent of τ and α. In our case, the following simple example of a function γ(τα) satisfying (8.8) is relevant: where √ 1 + τα denotes the principal branch of the square root function, compare [35]. Now, we have where γ(τα) satisfies (8.8). Then, there is a small number τ 0 > 0 such that for τ < τ 0 , (1) if γ(0) < 1, then there is a positive eigenvalue to (A 1); (2) if γ(0) > 1, then for any non-zero eigenvalue α of (8.9), we have Proof Lemma 13 follows from Theorem 19 by a regular perturbation argument. To make sure that the perturbation argument works, we have to show that if α R > −c (for some c > 0) and 0 6 τ < τ 0 (for some τ 0 > 0), where α = α R + √ −1α I , then |α| 6 C, where C is a generic constant which is independent of τ. In fact, multiplying (8.9) by the conjugatē φ of φ and integration by parts, we obtain that From the imaginary part of (8.10), we obtain that where C 1 is a positive constant (independent of τ). By assumption (8.8), |γ(τα)| 6 C and so |α I | 6 C. Taking the real part of (8.10) and noting that l.h.s. of (8.10) > C R |φ| 2 for some C ∈ R, we obtain that α R 6 C 2 , where C 2 is a positive constant (independent of τ > 0). Therefore, |α| is uniformly bounded and hence a regular perturbation argument gives the desired conclusion.
In conclusion, we have finished the study of the large eigenvalues (of order O(1)) and derived results on their stability properties.
It remains to study the small eigenvalues (of order o(1)) which will be done in the next section.

Stability proof II: characterization of small eigenvalues
Now, we study the eigenvalue problem (8.1), (8.2) with respect to small eigenvalues. Namely, we assume that λ ε → 0 as max ε √ D , D → 0. We will show that that the small eigenvalues are given by The matrix M(t 0 ) will be defined in (9.8) and given to leading order in (9.12). Before defining and computing the matrix, we have to make a few preparations.
where t ε = (t ε 1 , . . . , t ε N ) ∈ Ω η . After re-scaling, the eigenvalue problem (8.1), (8.2) becomes Here and in the rest of the proof for small eigenvalues, we set τ = 0. Since for small eigenvalues, we have τλ ε → 0 the proof and results extended to the case of a fixed constant τ > 0. Throughout this section, we denote By the implicit function theorem, there exists a (locally) unique solutionξ(t) = (ξ 1 (t), . . . ,ξ N (t)) of the equation Moreover,ξ(t) is C 1 for t ∈ Ω η . Note that we do not want to consider the solutionξ(t) = 0 since it does not correspond to a strictly positive solution.
We have the estimateŝ As a preparation, we first compute the derivatives ofξ(t). Now from (9.3), we calculate Here, F(t) is the vector field We compute Comparing with (7.5) and Proposition 11, we have . In addition, if M(t ε )) is positive definite, then we will show that all small eigenvalues have negative real part when 0 6 τ < τ 0 for some τ 0 > 0.
Next, we compute M(t) using (9.4). For |i − j| = 1, we compute in case and a similar result holds for This implies Therefore, using (7.6), (7.7) and the estimate (9.7), we have The matrix M(t ε ) will be the leading-order contribution to the small eigenvalues (compare Lemma 23 and the comments following it). Thus, we study the spectrum of the symmetric The corresponding eigenvectors are computed recursively from (9.17).
The matrix A has eigenvalue λ 1 = 0 with eigenvector v 1 = e. To compute the other eigenvalues and eigenvectors of A, we remark that this problem is equivalent to finding a suitable finite-difference approximationũ of the differential equation in the interval (0, 1) for uniform step-size h = 1 N . More precisely, we identify To determine the eigenvectors v i , we have to solve this finite-difference problem exactly. We assume that the solutions are given by polynomials of degree n (which will be shown later and n will be specified). Using Taylor expansion around x = x k−1/2 and the identities the finite-difference problem is equivalent to Substituting the ansatzũ a k x k into this equation, considering the coefficient of the power x k , k = 0, . . . , n, implies that (λ n − k(k + 1))a k + (k + 1) 2 a k+1 where for k = 0, we put (0 − 1)! = 1 in the second line of (9.17). For k = n, n = 0, 1, . . . , N − 1, this gives (λ n − n(n + 1))a n + (n + 1) 2 a n+1 = 0.
Case 2. n > N: Then, v n = 0 althoughũ ) 0. The resulting eigenfunctions for A are trivial and so in this case there are no new eigenpairs. Thus, we have found N eigenpairs with linearly independent eigenvectors.

Remark 17
The eigenvector v 0 with eigenvalue λ 0 = 0 corresponds to a rigid translation of all N spikes.
The leading eigenpair for mutual movement of spikes is (λ 1 , v 1 ).
The eigenvector for λ 1 = 2 can be computed as follows: The components of v 1,k are linearly increasing and have odd symmetry around the centre of the spike cluster which corresponds to k = N+1 2 or x = 1 2 .

Remark 18
The stability of the small eigenvalues follows from the results in [29] but the eigenvalues have not been determined explicitly.
The technical analysis for the small eigenvalues has been postponed to Appendix B.

Conclusion
We end this paper with a discussion of our results. We have considered a particular biological reaction-diffusion system with two small diffusivities, the Gierer-Meinhardt system with precursor. We have proved the existence and stability of cluster solutions which have three different length scales: a scale of order O(1) coming from the precursor inhomogeneity and two small scales which are of the same size as the square roots of the small diffusivities. In particular, the cluster solution can be stable for a suitable choice of parameter values. Such systems and their solutions play an important role in biological modelling to account for the bridging of length scales, e.g. between genetic, nuclear, intra-cellular, cellular and tissue levels. Our solutions incorporate and combine multiple scales in a robust and stable manner. A particular example of biological multi-scale patterns concerns the pattern formation of head (more precisely, hypostome), tentacles and foot in hydra. Meinhardt's model [17] correctly describes the following experimental observation: With tentacle-specific antibodies, Bode et al. [3] have shown that after head removal tentacle activation first reappears at the very tip of the gastric column. Then, this activation becomes shifted away from the tip to a new location, where the tentacles eventually appear. There are different lengthscales involved for this tentacle pattern: diameter of the gastric column, distance between tentacles and diameter of tentacles.
Let us describe the relation of this paper to [17] in more detail. The model in [17] can be explained in simplified form as follows: It consists of three activator-inhibitor systems, accounting for the formation of head, foot and tentacles, respectively. These subsystems are coupled by a joint source density. Further, there is direct interaction between the tentacle and head components to account for suppression of tentacle peaks at the site of head peaks. Altogether, the model is a seven-component reaction-diffusion system.
The main link to the results in this paper is to understand tentacle activation near a maximum of the source density. For this effect, the foot components can be neglected and so we are dealing with a five component system only. It is observed experimentally in [3] and computed numerically in [17] that near a sufficiently high local maximum of the source density tentacle peaks appear. Two different cases are studied: (i) if there is already a peak of head activator at this position, the tentacle peaks will appear at ring-shaped positions with the head activator peak in the centre of the ring, or (ii) if there is no previous head activator peak at this position, a tentacle activator peak will form, followed by a head activator peak which causes the tentacle peak to split into multiple peaks which are finally displaced to positions in ring-shaped positions with the activator peak in its centre.
For this effect to happen, it is assumed that the source density changes very slowly in time and acts on a rather long length scale. This corresponds to the precursor inhomogeneity in our model which is independent of time and has an O(1) length scale. The way the source density enters into the model is set up differently than in our paper resulting in a local maximum in [17] having a similar effect to a local minimum in our paper. We try to model some of these phenomena in a "minimal model" which consists of only two components corresponding to the two tentacle components in [17] coupled to a time-independent source density acting on an O(1) length scale.
The cluster pattern of spikes located in a sub-interval studied in our paper resembles the ring of tentacle peaks reduced to one dimension. (Work on the two-dimensional case is currently in progress.) It is interesting to note that our paper is successful in modelling isolated tentacles (without head formation) observed in some experimental situations as discussed in [17].
Comparing the two models in this paper and [17] leads to immediate possible extensions of the spike cluster analysis to models which are biologically more realistic by taking into account the following phenomena: (i) the effect of the head activator-inhibitor system could be added to show that the head activator peak pushes out the tentacle activator peaks; (ii) the tentacle activator peaks split easily due to saturation non-linearities in the tentacle subsystem (whereas in this paper, we do not consider the effect of saturation); (iii) add the foot activator-inhibitor system; (iv) replace the precursor inhomogeneity by a time-dependent source density which interacts with the other subsystems dynamically, e.g. it is enhanced by head activator, suppressed by foot activator, diffuses and possibly has its own predetermined inhomogeneity.
There are links of the model in [17] to other fields in biology such as the periodic spacing of secondary structures around a primary organizing region which is observed in the arrangement of leaves and flower elements in plants around the primary meristem [5].
The molecular basis underlying the model in [17] has recently been confirmed experimentally: After treatment of hydra with Alsterpaullone (which stabilizes β-catenin and thus increases the source density), it has been found that tentacle formation occurs over the whole body column [4]. Numerical computations have confirmed this behaviour [19]. This is in agreement with the pattern of multiple spikes covering the whole interval computed in our paper (see .
Systems of the type considered in this paper are a key to understanding the hierarchy of multi-stage biological processes such as in signalling pathways, where typically first largescale structures appear which induce patterns on successively smaller scales. The precursor can represent previous information from an earlier stage of development leading to the formation of fine structure at the present time. The multi-spike cluster in this paper is a typical small-scale pattern which is established near a pre-existing large-scale precursor inhomogeneity. where where u ∈ H 2 (R). Then, the conjugate operator of L under the scalar product in L 2 (R) is given by Then, we have the following result. To prove (A 6), we proceed in a similar way for L * . The lth equation of (A 4) is given as follows: Multiplying (A 9) by w and integrating, we obtain R w 2 Ψ l dy = 0.
Proof This result follows from the Fredholm Alternative and Lemma 20.
Finally, we study the eigenvalue problem for L: We have