Multi-spike Patterns in the Gierer–Meinhardt System with a Nonzero Activator Boundary Flux

Abstract

The structure, linear stability, and dynamics of localized solutions to singularly perturbed reaction–diffusion equations have been the focus of numerous rigorous, asymptotic, and numerical studies in the last few decades. However, with a few exceptions, these studies have often assumed homogeneous boundary conditions. Motivated by the recent focus on the analysis of bulk-surface coupled problems, we consider the effect of inhomogeneous Neumann boundary conditions for the activator in the singularly perturbed one-dimensional Gierer–Meinhardt reaction–diffusion system. We show that these boundary conditions necessitate the formation of spikes that concentrate in a boundary layer near the domain boundaries. Using the method of matched asymptotic expansions, we construct boundary layer spikes and derive a class of shifted nonlocal eigenvalue problems analogous to those studied in Maini et al. (Chaos 17(3):037106, 2007) for which we rigorously prove partial stability results. Moreover, by using a combination of asymptotic, rigorous, and numerical methods we investigate the structure and linear stability of selected one- and two-spike patterns. In particular, we find that inhomogeneous Neumann boundary conditions increase both the range of parameter values over which asymmetric two-spike patterns exist and are stable.

Appendices

Appendix A. Large $\lambda _I$ Asymptotics of ${\mathscr {F}}_{y_0}(i\lambda _I)$

Fig. 11

a Plot of ${\mathscr {F}}_{y_0}(0)$ versus the shift parameter $y_0$. b and c Real and imaginary parts of ${\mathscr {F}}_{y_0}(i\lambda _I)$ for select values of $y_0\ge 0$. The dashed lines indicate the $\lambda _I\gg 1$ asymptotics

In this appendix, we determine some key properties of ${\mathscr {F}}_{y_0}(\lambda )$ defined in (3.12). Recalling (3.13), in Fig. 11a we plot ${\mathscr {F}}_{y_0}(0)$ versus $y_0\ge 0$. Next, we calculate the limiting behaviour of ${\mathscr {F}}_{y_0}(i\lambda _I)$ as $\lambda _I\rightarrow \infty$. First, we let $({\mathscr {L}}_{y_0} - i\lambda _I)^{-1} w_c(y+y_0)^2 = \Phi _R + i\Phi _I$ where $\Phi _R$ and $\Phi _I$ solve

$$\begin{aligned} {\mathscr {L}}_{y_0}\Phi _R + \lambda _I\Phi _I = w_c(y+y_0)^2,\qquad {\mathscr {L}}_{y_0}\Phi _I - \lambda _I\Phi _R = 0, \end{aligned}$$

(A.1)

with the boundary conditions $\Phi _R'(0) = \Phi _I'(0)$ and $\Phi _R,\Phi _L\rightarrow 0$ as $y\rightarrow \infty$. Taking $\lambda _I\gg 1$ and assuming that $y={\mathcal {O}}(1)$, we obtain

$$\begin{aligned} \Phi _I(y) \sim \frac{1}{\lambda _I}w_c(y+y_0)^2,\qquad \Phi _R(y) \sim \frac{1}{\lambda _I}{\mathscr {L}}_{y_0}\Phi _I = \frac{1}{\lambda _I^2}\bigl ( 2w_c'(y+y_0)^2 + w_c(y+y_0)^2\bigr ). \end{aligned}$$

If $y_0>0$ then $\Phi _R'(0)=0$ and $\Phi _I'(0)=0$ are not satisfied and we must therefore consider the boundary layer at $y=0$. Setting $z = \lambda _I^{1/2} y$, we consider the inner expansion $\Phi _R \sim {\tilde{\Phi }}_R(z)$ and $\Phi _I \sim {\tilde{\Phi }}_I(z)$ where ${\tilde{\Phi }}_I$ satisfies

$$\begin{aligned} \frac{\text {d}^4{\tilde{\Phi }}_I}{\text {d}z^4} + {\tilde{\Phi }}_I = \frac{1}{\lambda _I} w_c(y_0)^2\quad z > 0;\qquad \frac{\text {d}{\tilde{\Phi }}_I}{\text {d}z} = \frac{d^3 {\tilde{\Phi }}_I}{dz^3} = 0,\quad z=0, \end{aligned}$$

and must be matched to the outer, $y={\mathcal {O}}(1)$, solution through the far-field behaviour

$$\begin{aligned} {\tilde{\Phi }}_I \sim \frac{1}{\lambda _I}w_c(y_0)^2,\quad \frac{\text {d}^2{\tilde{\Phi }}_I}{\text {d}z^2} \sim \frac{1}{\lambda _I^2}\bigl (2w_c'(y+y_0)^2 + w_c(y+y_0)^2\bigr ),\qquad z\rightarrow \infty . \end{aligned}$$

It is clear that the leading order solution is ${\tilde{\Phi }}_I(z) \sim \lambda _I^{-1} w_c(y_0)^2$. The constant behaviour of $\Phi _I$ at the boundary layer therefore does not contribute to the leading order behaviour of the integral

$$\begin{aligned} \int _0^\infty w_c(y+y_0)\Phi _I(y) {\text {d}}y \sim \lambda _I^{-1} \int _0^\infty w_c(y+y_0)^3 {\text {d}}y,\qquad \lambda _I\gg 1. \end{aligned}$$

Moreover, multiplying the right equation in (A.1) by $w_c(y+y_0)$ and integrating we calculate

$$\begin{aligned} \int _0^\infty w_c(y+y_0)\Phi _R(y) {\text {d}}y&= \frac{1}{\lambda _I}\biggl ( w_c'(y_0)\Phi _I(0) + \int _0^\infty \Phi _I{\mathscr {L}}_{y_0} w_c(y+y_0) {\text {d}}y \biggr )\\&\sim \frac{1}{\lambda _I^2}\biggl (w_c'(y_0)w_c(y_0)^2 + \int _0^\infty w_c(y+y_0)^4 {\text {d}}y\biggr ), \end{aligned}$$

for $\lambda _I\gg 1$ where we have used ${\mathscr {L}}_{y_0} w_c(y+y_0) = w_c(y+y_0)^2$. In summary, we have the large $\lambda _I$ asymptotics

$$\begin{aligned} {\mathscr {F}}_{y_0}(i\lambda _I) \sim \frac{1}{\lambda _I^2}\frac{w_c'(y_0)w_c(y_0)^2 + \int _0^\infty w_c(y+y_0)^4dy}{\int _0^\infty w_c(y+y_0)^2{\text {d}}y} + \frac{i}{\lambda _I}\frac{\int _0^\infty w_c(y+y_0)^3 {\text {d}}y}{\int _0^\infty w_c(y+y_0)^2{\text {d}}y},\qquad \lambda _I\gg 1. \end{aligned}$$

(A.2)

Note that the real part changes from positive to negative as $y_0$ exceeds $y_0\approx 1.0487$. In Fig. 11b, c we plot the real and imaginary parts of ${\mathscr {F}}_{y_0}(i\lambda _I)$, respectively, for select values of $y_0$. In addition, we have included the large $\lambda _I$ asymptotics which indicate close agreement for moderately large values of $\lambda _I$.

Appendix B. Numerical Support for Stability Conjecture

In this appendix, we provide numerical support for the conjecture that the shifted NLEP (4.1) has a stable spectrum when $\mu >\mu _c(y_0)$ by numerically calculating the dominant eigenvalue of the NLEP for $0\le y_{0} \le 1.5$ and $0\le \mu \le 10$. The numerical calculation of the spectrum was performed by truncating the domain $0<y<\infty$ to $0<y<20$ and discretizing it with 600 uniformly distributed points. Then, we used a finite difference approximation for the second derivatives and a trapezoidal rule discretization for the integral term to approximate the NLEP (4.1) with a discrete matrix eigenvalue problem. We then numerically calculated the dominant eigenvalue of matrix by using the eigenfunction in the Python scipy.linalg library for our numerical computation of the dominant eigenvalue. In Fig. 12a we plot $\text {Re}\lambda _0$ versus $y_0$ and $\mu$. We observe the real part of the dominant eigenvalue is negative when $\mu$ exceeds the threshold $\mu _c(y_0)$. Additionally, in Fig. 12b we plot $\Lambda _0-\text {Re}(\lambda _0)$ for the same range of $y_0$ and $\mu$ values. We observe that this difference is non-negative which suggest that $\text {Re}\lambda _0\le \Lambda _0$.

Fig. 12

a Plot of the real part of the dominant eigenvalue of the shifted NLEP (4.1) versus shift parameter $y_0$ and multiplier $\mu$. The dotted red line corresponds to the critical threshold $\mu _c$ defined in (4.14) and the solid dark line is the zero contour of $\text {Re}\lambda _0$. b Plot of the difference between dominant eigenvalues of ${\mathscr {L}}_{y_0}$ and the NLEP (4.1)

Appendix C. Stability of Asymmetric Two-Boundary-Spike Pattern when $A=0$

Previous results on the stability of asymmetric two spike equilibria of (1.2) when $A=0$ have focused exclusively on interior multi-spike solutions (Ward and Wei 2002). To compare the $A=B=0$ theory with our results obtained in Examples 2–4, we include here a summary of the stability of an asymmetric two-spike solution where one spike concentrates at $x=0$ and the other concentrates at either $x=1$ or in the interior $0<x<1$. In both cases the NLEP (3.6a) with $A=0$ can be written as

$$\begin{aligned} {\mathscr {L}}_0\pmb {\phi } - 2 w_c(y)^2 \frac{\int _0^\infty w_c(y) {\mathscr {E}}\pmb {\phi }{\text {d}}y}{\int _0^\infty w_c(y)^2 {\text {d}}y} = \lambda \pmb {\phi }, \end{aligned}$$

(C.1)

where for two boundary spikes we let

Fig. 13

a Multiplier $\mu$ in the NLEP for a two-spike solution when $A=B=0$ and where the two spikes concentrate either at the boundaries or else one spike concentrates at the boundary and the other in the interior. b Relationship between l and D in the gluing method used to construct a two-spike solution when $A=B=0$ and where either both spikes concentrate at the boundaries, or one spike concentrates at the boundary and the other in the interior

$$\begin{aligned} {\mathscr {E}} = {\mathscr {E}}_{bb} \equiv \begin{pmatrix} \coth \omega _0\tanh \omega _0l &{} \text {csch}\omega _0 \tanh \omega _0(1-l) \\ \text {csch}\omega _0\tanh \omega _0l &{} \coth \omega _0\tanh \omega _0(1-l) \end{pmatrix}, \end{aligned}$$

(C.2)

and in the case of one boundary and one interior spike we let

$$\begin{aligned} {\mathscr {E}} = {\mathscr {E}}_{bi} \equiv \begin{pmatrix} \coth \omega _0\tanh \omega _0l &{} 2\text {csch}\omega _0\sinh \omega _0\tfrac{1-l}{2} \\ \text {csch}\omega _0\tanh \omega _0l\cosh \omega _0\tfrac{1-l}{2} &{} 2\text {csch}\omega _0\cosh \omega _0\tfrac{1+l}{2}\sinh \omega _0\tfrac{1-l}{2} \end{pmatrix}. \end{aligned}$$

(C.3)

It is then straightforward to verify that $\sigma =1$ is an eigenvalue of both matrices ${\mathscr {E}}_{bb}$ and ${\mathscr {E}}_{bi}$. The remaining eigenvalue in each case is then given by the determinant. By diagonalizing ${\mathscr {E}}$ the NLEP (C.1) can therefore be written as (4.1) with $\mu =2$ as well as with $\mu =2\det {\mathscr {E}}_{bb}$ and $\mu = 2\det {\mathscr {E}}_{bi}$ for the boundary–boundary and boundary–interior cases, respectively. Since the $A=B=0$ stability theory implies that the NLEP (4.1) is stable if and only if $\mu >1$ (Wei 1999) we immediately deduce that the $\mu =2$ modes are stable in both the boundary–boundary and boundary–interior cases. To determine the stability of the remaining modes, we explicitly calculate

$$\begin{aligned} \det {\mathscr {E}}_{bb} = \tanh \omega _0 l\tanh \omega _0(1-l),\qquad \det {\mathscr {E}}_{bi} = 2\frac{\tanh \omega _0l\sinh \omega _0\tfrac{1-l}{2}\sinh \omega _0\tfrac{1+l}{2}}{\sinh \omega _0}. \end{aligned}$$

(C.4)

Finally, for the boundary–boundary and boundary–interior cases we solve (5.7) and (5.20) for $D=D(l)$, respectively, and the resulting values of $2\det {\mathscr {E}}_{bb}$ and $2\det {\mathscr {E}}_{bi}$ for $0<l<1$ are shown in Fig. 13a. In particular, the asymmetric pattern with two boundary spikes is always linearly unstable, while the pattern with one boundary and one interior spike has a region of stability (with respect to the ${\mathcal {O}}(1)$ eigenvalues). In Fig. 13b we plot $l=l(D)$ (cf. Figs. 7 and 10) for both two-spike configurations, indicating where the pattern is stable (solid line) and unstable (dashed line).