Dynamics of shadow system of a singular Gierer-Meinhardt system on an evolving domain

The main purpose of the current paper is to contribute towards the comprehension of the dynamics of the shadow system of a singular Gierer-Meinhardt model on an isotropically evolving domain. In the case where the inhibitor's response to the activator's growth is rather weak, then the shadow system of the Gierer-Meinhardt model is reduced to a single though non-local equation whose dynamics is thoroughly investigated throughout the manuscript. The main focus is on the derivation of blow-up results for this non-local equation, which can be interpreted as instability patterns of the shadow system. In particular, a {\it diffusion-driven instability (DDI)}, or {\it Turing instability}, in the neighbourhood of a constant stationary solution, which then is destabilised via diffusion-driven blow-up, is observed. The latter indicates the formation of some unstable patterns, whilst some stability results of global-in-time solutions towards non-constant steady states guarantee the occurrence of some stable patterns. Most of the derived results are confirmed numerically and also compared with the ones in the case of a stationary domain.


Introduction
The purpose of the current work is to study an activator-inhibitor system, introduced by Gierer and Meinhard in 1972 [5] to describe the phenomenon of morphogenesis in hydra, on an evolving domain. Assume that u(x, t) stands for the concentration of the activator, at a spatial point x ∈ Ω t ⊂ R N , N = 1, 2, 3, at time t ∈ [0, T ], T > 0, which enhances its own production and that of the inhibitor. On the other hand, let v(x, t) represents the concentration of the inhibitor, which suppresses its own production as well as that of the activator. Hence, the interaction between u and v can be described by the following non-dimensionalised system [5] u t + ∇ · ( − → α u) = D 1 ∆u − u + u p v q , x ∈ Ω t , t ∈ (0, T ), (1.1) where ν is the unit normal vector on ∂Ω t , whereas − → α ∈ R N stands for the convection velocity which is induced by the material deformation due to the evolution of the domain. Moreover, D 1 , D 2 are the diffusion coefficients of the activator and inhibitor respectively; τ represents the response of the inhibitor to the activator's growth. Moreover, the exponents satisfying the conditions: p > 1, q, r, > 0, and s > −1, measure the interactions between morphogens. The dynamics of system (1.1)-(1.4) can be characterised by two values: the net self-activation index π = (p − 1)/r and the net cross-inhibition index γ = q/(s + 1). Index π correlates the strength of self-activation of the activator with the cross-activation of the inhibitor. Thus, if π is large, then the net growth of the activator is large no matter the growth of the inhibitor. The parameter γ measures how strongly the inhibitor suppresses the production of the activator and that of itself. If γ is large then the production of the activator is strongly suppressed by the inhibitor. Finally, the parameter τ quantifies the inhibitor's response against the activator's growth. Guided by biological interpretation as well as by mathematical reasons, we assume that the parameters p, q, r, s satisfy the condition p − rγ < 1, (1.5) which in the literature is known as the Turing condition since it guarantees the occurrence of Turing patterns for the system (1.1)-(1.4) on a stationary domain [14]. For analytical purposes, in the current work we will only consider the case of an isotropic flow on an evolving domain, and thus we have for any x ∈ Ω t : with ρ(t) being C 1 −function with ρ(0) = 1. In the case of a growing domain we havė ρ(t) = dρ dt > 0, whilst when the domain shrinks or for domain contractionρ(t) = dρ dt < 0. Furthermore, the following equality holds (1.7) 2 Settingû(ξ, t) = u(ρ(t)ξ, t),v(ξ, t) = v(ρ(t)ξ, t), and then using the chain rule as well as (1.6) and (1.7), see also [12], we obtain: whilst similar relations hold for v as well. Therefore (1.1)-(1.4) is reduced to following system on a reference stationary domain Ω 0 where ∆ ξ represents the Laplacian on the reference static domain Ω 0 . Henceforth, without any loss of generality we will omit the index ξ from the Laplacian. Defining a new time scale [10], 12) and settingũ(ξ, σ) =û(ξ, t),ṽ(ξ, σ) =v(ξ, t), then system (1.8)-(1.11) can be written as where ρ(t) = φ(σ), and thusρ(t) =φ (σ) φ 2 (σ) , and Σ = σ(T ). Now if D 1 D 2 , i.e. when the inhibitor diffuses much faster than the activator, then system (1.13)-(1.16) can be fairly approximated by an ODE-PDE system with a nonlocal reaction term. We will denote the new approximation by shadow system as coined in [9]. Below we provide a rather rough derivation of the shadow system, while for a more rigorous approach one can appeal to the arguments in [1]. Indeed, dividing (1.14) by D 2 and taking D 2 → +∞, see also [15], then it follows thatṽ solves ∆ ξṽ = 0, ξ ∈ Ω 0 , ∂ṽ ∂ν = 0, ξ ∈ ∂Ω 0 , for any fixed σ ∈ (0, Σ). Due to the imposed Neumann boundary condition thenṽ is a spatial homogeneous (independent of ξ) solution, i.e.ṽ(ξ, σ) = η(σ) and thus (1.14) can be written as . (1.18) Averaging (1.17) over Ω 0 we finally infer that the pair (ũ, η) satisfies the shadow system where − Ω 0ũ r dξ = 1 |Ω 0 | Ω 0ũ r dξ. In the limit case τ → 0, i.e. when the inhibitor's response to the growth of the activator is quite small, then the shadow system is reduced to a single, though, non-local equation. Indeed, when τ = 0, (1.20) entails that η(σ) = recalling γ = q s+1 and Recovering the t variable entails that the following partial differential equation holdŝ The main aim of the current work is to investigate the long-time dynamics of the non-local problem (1.23)-(1.25) and then check whether it resembles that of the reactiondiffusion system (1.13)-(1.16). Biologically speaking we will investigate whether, under the fact that the inhibitor's response to the growth of the activator is quite small and when it also diffuses much faster than the activator, is it necessary to study the dynamics of both reactants or it is sufficient to study only the activator's dynamics. From here onwards, we take D 1 = 1, revert to the initial variables x, u instead of ξ, u and we drop the index ξ from the Laplacian ∆ without any loss of generality. Hence, we will focus our study on the following single partial differential equation The layout of the current work is as follows. Section 2 deals with the derivation and proofs of various blow-up results, induced by the non-local reaction term (ODE blow-up results), together with some global-time existence results for problem (1.30)-(1.25). Following the approach developed in [7,8], in Section 3 we present and prove a Turing instability result associated with (1.

ODE Blow-up and Global Existence
The current section is devoted to the presentation of some blow-up results for problem (1.30)-(1.32), i.e. blow-up results induced by the kinetic (non-local) term in (1.30). Besides, some global-in-time existence results for problem (1.30)-(1.32) are also presented. Throughout the manuscript we use the notation C and c to denote positive constants with big and small values respectively. Our first observation is that the concentration of the activator cannot extinct in finite time. Indeed, the following proposition holds.
then for each Σ > 0 there exists C Σ > 0 such that for the solution u(x, σ) of (1.30)-(1.32) the following inequality holds Proof. Owing to the maximum principle and by using (2.1) we derive that u = u(x, σ) > 0. By virtue of the comparison principle, we also deduce that u(x, σ) ≥ũ(σ), whereũ =ũ(σ) is the solution to dũ Remark 2.1. It is easily checked that condition (2.1) is satisfied for any decreasing function φ(σ) satisfying since then by virtue of (1.18) A key estimate for obtaining some blow-up results presented throughout is the following proposition.
Remark 2.2. Note that Proposition 2.2 guarantees that the non-local term of problem (1.30)-(1.32) stays away from zero and hence its solution u is bounded away from zero as well. In fact, inequality (2.7) implies −

14)
follows by Jensen's inequality, [3], taking δ ≤ r, where again c is independent of time t. The latter estimate rules out the possibility of (finite or infinite time) quenching, i.e. lim σ→Σ ||u(·, σ)|| ∞ = 0, cannot happen, and thus extinction of the activator in the long run is not possible.

Remark 2.3.
In case Φ(σ) is not bounded from above, as it happens for ρ(t) = e βt , β > 0, , then both of the estimates (2.7) and (2.14) still hold true, however the involved constants depend on time σ and thus (finite or infinite time) quenching cannot be ruled out.
Next we present our first ODE-type blow-up result for problem (1.30)-(1.32) when an anti-Turing condition, the reverse of (1.5), is satisfied.

19)
and In that case and again finite-time blow-up occurs at For a stationary domain, i.e. when ρ(t) = φ(σ) = 1, we have Φ(σ) = Ψ(σ) = 1 and thus finite-time blow-up occurs at Remark 2.5. When the domain evolves logistically, a feasible choice in the context of biology [16], i.e. when ρ(t) = e βt 1+ 1 m (e βt −1) , for m = 1, means that (1.12) cannot be solved for t and it is more convenient to deal with problem (1.27)-(1.29) instead. Then following the same approach as in Theorem 2.1 we show that the solution of (1.27)-(1.29) exhibits finite-time blow-up under the same conditions for parameters p, γ, r provided that the initial condition satisfies where now the quantity L(t) = 1 + Remark 2.6. Assume now that then G(Σ) > 0 and since G(σ) is strictly decreasing we get that G(σ) > 0 for any 0 < σ < Σ which implies that F (σ) never blows up. Therefore, since F (σ) ≤ū(σ), there is still a possibility thatū(σ) does not blow up either, however we cannot be sure and it remains to be verified numerically, see Section 4.
Next, we investigate the dynamics of some L -norms ||u(·, σ)|| , which identify some invariant regions in the phase space. We first define Our first result in this direction provides some conditions under which a finite-time blowup takes place, when an anti-Turing condition is in place, and is stated as follows.

Turing Instability and Pattern Formation
In the current section, we present a Turing-instability result for problem (1.30)-(1.32), restricting ourselves to the radial case Ω 0 = B 1 (0) := {x ∈ R N | |x| < 1}. Then the solution of u (1.30)-(1.32) is radially symmetric, i.e. u(x, t) = u(R, t) for R = |x|, and thus it satisfies the following where ∆ R u := u RR + N −1 R u R . It can be seen, see also [7,8], that under the Turing condition (1.5), the spatial homogeneous solutions of (1.30)-(1.32), i.e. the solution of the problem du dσ = −Φ(σ)u + Ψ(σ)u p−rγ , u| σ=0 =ū 0 > 0, never exhibits blow-up, as long as Φ(σ), Ψ(σ) are both bounded, since the nonlinearity f (u) = u p−rγ is sub-linear. On the other hand, considering spatial inhomogeneous solutions of (1.30)-(1.32) we will show below, see Theorem 3.1, that a diffusion-driven instability phenomenon occurs when spiky type of initial conditions are considered. Indeed, next we consider an initial datum of the form, [6], with 0 < λ 1 and where a = 2 p−1 and 0 < δ < 1. It is easily seen, [7,8], that for ψ δ the following lemma holds. Furthermore, it follows that (3.10) and the initial data u 0 defined by (3.4) also satisfies the following lemma, see also [7,8], Hereafter, we fix 0 < λ ≤ λ 0 = λ 0 (d) so that (3.11) is satisfied. Given 0 < δ < 1, let Σ δ > 0 be the maximal existence time of the solution to (3.1)-(3.3) with initial data of the form (3.4)-(3.5). Next, we introduce the new variable z = e σ Φ(s) ds u, so that the linear dissipative term −Φ(σ)u is eliminated and then z satisfies (3.13) z(R, 0) = u 0 (R), 0 < R < 1, (3.14) where It is clear that as long as Φ(σ) is bounded then u blows-up in finite time if and only if z does so. Assuming now that both Φ(σ) and Ψ(σ) are positive and bounded, which is the case for the evolution provided by ψ(σ) satisfying (2.3) or for an exponential shrinking domain as indicated in Remarks 2.1 and 2.4, then by virtue of (2.14) we have thus averaging of (3.12) entails Henceforth, the positivity and the boundedness of Φ(σ), Ψ(σ) as well as the Turing condition (1.5) are imposed. Moving towards the proof of Theorem 3.1 we first need to establish some auxiliary results. Next, we provide a useful estimate of z that will be frequently used throughout the text.
for any 0 < δ < 1 and some positive constant c.
Next note that when p > N N −2 , there is 1 < q < p such that N > 2p q−1 and thus the following quantities are finite due to (3.9). An essential ingredient for the proof of Theorem 3.1 is the following key estimate of the L p −norm of z in terms of A 1 and A 2 .
The proof of Proposition 3.1 requires some further auxiliary results provided below. Let us define 0 < σ 0 (δ) < Σ δ to be the maximal time for which inequality (3.25) is valid in 0 < σ < σ 0 (δ), then we have We only regard the case σ 0 (δ) ≤ 1, since otherwise there is nothing to prove. Then the following lemma holds true.
On the other hand, for δ ≤ R ≤ 1 and by using (3.7) for m = 1 we take for 0 < σ < min{σ 0 , Σ δ }, and therefore, (3.45) Note that for 0 < δ 1, then the right-hand side on (3.45) is less than σ 0 , so Σ δ < Remark 3.3. Notably, by (3.45) we conclude that Σ δ → 0 as δ → 0, i.e. the more spiky initial data we consider then the faster the diffusion-driven blow-up for z and thus for u take place.
A diffusion-driven instability (Turing instability) phenomenon, as was first indicated in the seminal paper [21], is often followed by pattern formation. A similar situation is observed as a consequence of the driven-diffusion finite-time blow-up provided by Theorem 3.1, and it is described below. The blow-up rate of the solution u of (3.1)-(3.3) and the blow-up pattern (profile) identifying the formed pattern are given.
Assume that both Φ(σ) and Ψ(σ) are positive and bounded. Then the blow-up rate of the solution of (3.1)-(3.3) can be characterised as follows

46)
where Σ max stands for the blow-up time.
Remark 3.4. We first perceive that (3.48) provides a rough form of the blow-up pattern for z and thus for u as well. Additionally, owing to (3.47) the non-local problem (3.12)-(3.14) can be treated as a local one for which the more accurate asymptotic blow-up profile, [13], is known and is given by lim σ→Σmax z(|x|, σ) ∼ C | log |x|| |x| 2 for |x| 1, and C > 0.
Using again the relation between z and u we end up with a similar asymptotic blow-up profile for the diffusion-driven-induced blow-up solution u of problem (3.1)- (3.3). This blow-up profile actually determines the form of the developed patterns which are induced as a result of the diffusion-driven instability and it is numerically investigated in the next section.

Numerical Experiments
To illustrate some of the theoretical results of the previous sections we perform a series of numerical experiments for which we solve the involved PDE problems using the Finite Element Method [18], using piecewise linear basis functions and implemented using the adaptive finite-element toolbox ALBERTA [19]. In all our simulations (unless stated otherwise) the domain was triangulated using 8192 elements, the discretisation in time was done using the forward Euler method taking 5 × 10 −4 as time-step and the resulting linear system solved using Generalized Minimal Residual iterative solver [20].  29) on Ω 0 = [0, 1] 2 . The initial condition is chosen u 0 (x, 0) = cos(πy) + 2. As for the domain evolution we consider four different cases: • ρ(t) = e βt (exponentially growing domain); • ρ(t) = e −βt (exponentially decaying domain); • ρ(t) = e βt 1+ 1 m (e βt −1) (logistically growing domain); • ρ(t) = 1 (static domain). We summarise all parameters used in Table 4.1. In Fig. 4.1, we demonstrate the ||u(x, t)|| ∞ for each of the domain evolutions, so we can monitor their respective blow-up times. If we denote by Σ g , Σ s , Σ lg and Σ 1 the blow-up times for the case of exponentially D 1 p q r s β m 1 3 2 1 2 0.1 1.5 Table 1. Set of parameters used in Experiment 1. growing and decaying, the logistically growing domains and the static domain, respectively, we observe from Fig. 4.1 that we have the following relation Σ g > Σ lg > Σ 1 > Σ s , which is in agreement with the mathematical intuition.
We now take the same initial condition, u 0 and the same initial domain which we assume is evolving exponentially and consider parameters D 1 = 1, p = 1.4, q = 1, r = 1 and s = 2 for which inequality (2.23) of Remark 2.6 holds. As we can see in Fig. 4.1, we have an example of a solutionū that does not blow up, as already conjectured in the aforementioned remark. Notably, this numerical experiment predicts a very interesting phenomenon both mathematically and biologically. It predicts the infinite-time quenching of the solution of problem (1.27)-(1.29) and thus the extinction of the activator in the long run, see also Remark 2.3. Note also that this result is not in contradiction with Proposition 2.2, where infinite-time quenching is ruled out since condition (2.1) is not satisfied for an exponentially growing domain where Φ(σ) is an unbounded function as indicated in Remark 2.4.

Experiment 2.
This experiment is meant to illustrate Theorem 2.3 and we take the same initial data u 0 = cos(πy) + 2 and take Ω 0 as the unit square when numerically solving equations (1.27)-(1.29). As for domain evolution we consider ρ(t) = e βt , with β = 0.1. To proceed, we consider two sets of parameters, one for which assumptions of Theorem 2.3 are satisfied and another for which those assumptions are not fulfilled. See Table 4 Results shown in Fig. 4.2 are in agreement with theoretical predictions of Theorem 2.3 since the solutions exists for all times when the assumption of the theorem are met (Fig. 4.2(a)), otherwise a finite-time blow-up is exhibited to occur (Fig. 4.2(b)).

Experiment 3.
In this experiment we intend to illustrate Theorem 3.1 so we numerically solve (1.27)-(1.29) in R 3 , taking Ω 0 as the unit sphere and initial condition u 0 given by (3.4), considering δ = 0.8 and λ = 0.1. As for other parameters we choose D 1 = 1, p = 4, q = 4, r = 2 and s = 1, which satisfy the conditions of the theorem. In Fig. 4.3 we display the L ∞ −norm of the solution u for three types of evolution laws implemented, namely: exponential decay, logistic decay and no evolution. For the exponential and logistic decay we select the same set of parameters as used in Experiment 1. As we can observe, for all the cases the solution blows up, as theoretically predicted by 25 Figure 3. The plot of ||u(x, t)|| ∞ , where u(x, t) is the numerical solution of (1.27)-(1.29). Initial condition is u 0 = cos(πy) + 2 and Ω 0 is the unit square evolving according to exponential growth (β = 0.1). (Colour version online). Theorem 3.1. Again the blow-up times have the order Σ 1 > Σ ls > Σ s , where now Σ ls stands for the blow-up time for the logistic decay evolution, beeing in agreement with the mathematical intuition. In Fig. 5(a) and Fig. 5(b) we compare the initial solution with the solution at t = 0.03 respectively, for the logistic decay, close to the blow-up time t = 0.03, by looking at a cross section of the three-dimensional unit sphere Ω 0 . Besides, in Fig. 5(c) and Fig. 5(d) again the solution at section cross of Ω 0 is depicted but now for the stationary and exponential decaying case respectively. Through this experiment we can observe the formation of blow-up (Turing-instability) patterns around the origin R = 0. We actually conclude that the evolution of the domain has no impact on the form of blow-up patterns, however it definitely affects the spreading of Turing-instability patterns as it is obvious from Fig. 5(b), Fig. 5(c) and Fig. 5(d). The obtained results are displayed in Fig. 4.4 and they actually demonstrate that reaction-diffusion system (1.13)-(1.14) and non-local problem (1.27)-(1.29) share the same dynamics. In particular the solutions of both problems exhibit blow-up which takes place in finite time. The latter, biologically speaking, means that in the examined case we just need to monitor only the dynamics of the activator, whose dynamics governed by nonlocal problem (1.27)-(1.29). Then we can get an insight regarding the interaction between both of the chemical reactants (activator and inhibitor) provided by reaction-diffusion system (1.13)-(1.14).