Construction and asymptotic behavior of multi-peak solutions to the Gierer-Meinhardt system with saturation

$A_t=\varepsilon^2 \Delta A-A+\frac{A^2}{H(1+kA^2)},\ A>0,\ $ in $\Omega\times (0,\infty), $

$\tau H_t=D\Delta H-H+A^2,\ H>0,\ $ in $\Omega \times (0,\infty),$

$\frac{\partial A}{\partial \nu}=\frac{\partial H}{\partial \nu}=0,\ $ on $\partial \Omega\times (0,\infty),$

where $\varepsilon>0$, $\tau \geq 0$, $k>0$. The unknowns $A=A(x,t)$, $H=H(x,t)$ represent the concentrations of the activator and the inhibitor at a point $x\in \Omega \subset R^N$ and at a time $t>0$. Here $\Delta$ := $\sum_{j=1}^N\frac{\partial^2}{\partial x^2_j}$ is the Laplace operator in $R^N$, $\Omega$ is a bounded smooth domain in $R^N$, and $\nu=\nu(x)$ is the outer unit normal at $x\in \partial \Omega$. When $\Omega$ is an $x_N$-axially symmetric domain and $2\leq N\leq 5$, for sufficiently small $\varepsilon>0$ and sufficiently large $D>0$ we construct a multi-peak stationary solution peaked at arbitrarily chosen intersections of $x^N$-axis and $\partial \Omega$, under the condition that $4k\varepsilon^{-2N}|\Omega|^2$ converges to some $k_0\in[0,\infty)$ as $\varepsilon\to 0$.