DELAYED REACTION KINETICS AND THE STABILITY OF SPIKES IN THE GIERER–MEINHARDT MODEL

. A linear stability analysis of localized spike solutions to the singularly perturbed two- component Gierer–Meinhardt (GM) reaction-diﬀusion (RD) system with a ﬁxed time delay T in the nonlinear reaction kinetics is performed. Our analysis of this model is motivated by the computational study of Lee, Gaﬀney, and Monk [ Bull. Math. Bio. , 72 (2010), pp. 2139–2160] on the eﬀect of gene expression time delays on spatial patterning for both the GM model and some related RD models. It is shown that the linear stability properties of such localized spike solutions are characterized by the discrete spectra of certain nonlocal eigenvalue problems (NLEP). Phase diagrams consisting of regions in parameter space where the steady-state spike solution is linearly stable are determined for various limiting forms of the GM model in both 1-dimensional and 2-dimensional domains. On the boundary of the region of stability, the spike solution is found to undergo a Hopf bifurcation. For a special range of exponents in the nonlinearities for the 1-dimensional GM model, and assuming that the time delay only occurs in the inhibitor kinetics, this Hopf bifurcation boundary is readily determined analytically. For this special range of exponents, the challenging problem of locating the discrete spectrum of the NLEP is reduced to the much simpler problem of locating the roots to a simple transcendental equation in the eigenvalue parameter. By using a hybrid analytical- numerical method, based on a parametrization of the NLEP, it is shown that qualitatively similar phase diagrams occur for general GM exponent sets and for the more biologically relevant case where the time delay occurs in both the activator and inhibitor kinetics. Overall, our results show that there is a critical value T (cid:63) of the delay for which the spike solution is unconditionally unstable for T > T ∗ , and that the parameter region where linear stability is assured is, in general, rather limited. A comparison of the theory with full numerical results computed from the RD system with delayed reaction kinetics for a particular parameter set suggests that the Hopf bifurcation can be subcritical, leading to a global breakdown of a robust spatial patterning mechanism.

1. Introduction. In [21], Turing proposed that localized peaks in the concentration of a chemical substance, known as a morphogen, could be responsible for the process of morphogenesis, which describes the development of a complex organism from a single cell. Through the use of a linearized analysis, he showed how stable spatially complex patterns can develop from small perturbations of spatially homogeneous initial data for a coupled two-component system of reaction-diffusion (RD) equations. Although there is now a vast literature on the study of Turing pattern formation for specific RD systems proposed to model various morphogenetic processes (cf. [5] and the references therein), these previous studies have typically neglected any time delays in the reaction kinetics owing to the time needed for gene expression. More specifically, as discussed in [5], there may exist a time delay between the initiation of protein signal transduction, due to ligand-receptor binding, and the time at which genes are ultimately produced.
In [5], [14], [15], [16] (see also the survey [18]), the effect of a fixed time delay in the reaction kinetics for some Turing pattern formation systems, by modeling various hypothetical subcellular gene expression dynamical processes, was studied computationally and through a Turing-type linear stability analysis on both fixed and slowly growing domains. Each of the two-component RD models studied in [5], [14], [15], and [16] is characterized by a short-range activation and a long-range inhibition, with the Gierer-Meihnardt (GM) RD model being one such prototypical system [6]. In [14], a GM model with a time-delayed reaction kinetics, modeling a signal transduction process involving reversible binding at the cell surface, was studied numerically (this is Model I in [14]). This computational study [14] showed that temporal oscillations in the spatial patterning occur as the time delay increases, and that these oscillations become large and uncontrolled as the delay increases further. The main goal of this paper is to provide a theoretical framework to predict parameter ranges where stable spatial patterning exists for this GM model of [14] with delayed reaction kinetics.
In one spatial dimension, the dimensionless GM RD model [6] allowing for delayed reaction kinetics (Model I of [14]), and posed on |x| ≤ 1 with u x (±1, t) = v x (±1, t) = 0, is 2 1 is the activator diffusivity, D > 0 is the inhibitor diffusivity, and τ > 0 is a reaction-time constant. We assume the usual condition (cf. [11]) on the exponents (p, q, m, s) that In the semistrong regime where 1 and D = O (1), and in the absence of delayed kinetics, there is a large literature on the stability of localized one-dimensional (1-D) spike solutions to the GM model (1.1) (cf. [25], [4], [10], [11], [24], [23], [22]). To analyze the linear stability of a steady-state spike solution to O(1) time-scale instabilities, the main technical challenge is that one must rigorously analyze the discrete spectrum of a nonlocal eigenvalue problem (NLEP). Although there are many rigorous results on the spectrum of the NLEP for various ranges of the exponent set (p, q, m, s) (see the survey [26]), the theory is still incomplete. More recently, in [20], it was shown that when p = 2m − 3 and m > 2, the study of the spectrum of the NLEP for the 1-D undelayed GM model (1.1) can be reduced to the analysis of a simple transcendental equation in the eigenvalue parameter. For this parameter range, p = 2m − 3 with m > 2, where we refer to the NLEP as being "explicitly solvable," detailed results for the linear stability of 1-D spike solutions are readily obtained. The specific case p = m = 3 also arises in the study of the stability of hot-spot patterns for an RD model of urban crime (cf. [13]).
Motivated by the previous computational studies (cf. [5], [14], [15], [16]) of spatial patterning for the GM and related models with delayed reaction kinetics, the main goal of this paper is to analyze the linear stability of steady-state spike solutions for (1.1) and its two-dimensional (2-D) counterpart to O(1) time scale instabilities when the reaction kinetics have a time delay T . Any such instability is an instability of the amplitude of the spike, and is associated with unstable O(1) eigenvalues in an NLEP. For various limiting forms of (1.1), as we discuss below, our main focus is to Downloaded 10/23/17 to 163.1.81.103. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php Figure 1, we plot the Hopf bifurcation boundary and the region of linear stability for the 1-D shadow problem for the prototypical GM exponent set (p, q, m, s) = (2, 1, 2, 0). We conclude that when T exceeds some threshold value the spike solution is always unstable.
In section 3 we analyze the linear stability of a one-spike steady state for the 1-D infinite-line problem for (1.1) assuming that the delay only occurs in the inhibitor kinetics. For this case, we set D = 1 and consider (1.1) on −∞ < x < ∞. For the parameter range p = 2m − 3 with m > 2 where the NLEP is explicitly solvable, in section 3.1 we analytically determine the Hopf bifurcation boundary in the τ versus T plane. For an arbitrary GM exponent set, in section 3.2 we formulate and implement a simple numerical method to readily compute the Hopf bifurcation boundary from a parametrization of the NLEP. For the prototypical GM exponent set (p, q, m, s) = (2, 1, 2, 0) this boundary and the region of linear stability is shown in the right panel of In section 4 we study the linear stability of an M -spot solution with M ≥ 2, for the GM model with prototypical parameter set (p, q, m, s) = (2, 1, 2, 0) in a 2-D bounded domain Ω with delayed inhibitor kinetics. In our analysis we consider the weakcoupling regime of [27], where the inhibitor diffusivity D satisfies D = O(− log ). In the absence of delay, it is well known (cf. [27]) that the multispot pattern can be destabilized by either a synchronous perturbation of the spot amplitudes or by M − 1 possible asynchronous perturbations of the spot amplitudes. By using a hybrid analytical-numerical approach to analyze the corresponding NLEP for either the synchronous or asynchronous modes, we determine phase diagrams in the threedimensional parameter space τ , T , and µ = 2πM D 0 /|Ω|, where the multispot pattern is linearly stable. Here D = D 0 /ν, ν ≡ −1/ log , and |Ω| is the area of Ω.
In section 5 we consider various limiting forms of the 1-D and 2-D GM model where we now assume that the reaction kinetics have a time delay in both the activator and the inhibitor. We first formulate a generalized NLEP to encompass the infiniteline problem, the shadow problem, and the multispot problem in two dimensions. Although it is analytically intractable to determine Hopf bifurcation boundaries from this NLEP, we provide a rigorous argument, based on continuity of paths of spectra, to prove that there is a critical value of the delay T at which a Hopf bifurcation occurs. Our numerical computations of the Hopf bifurcation boundary for the 1-D problem show that the phase diagram in the τ versus T plane is qualitatively very similar to the phase diagram that occurs when there is a time delay only in the inhibitor kinetics.
Finally, in section 6 we briefly summarize our main results and suggest a few open problems for further investigation.
2. Delay effects: The limiting shadow problem. In this section, we will analyze the shadow limit D → +∞ for the GM model (1.1). In the limit D → ∞, (1.1) reduces to the so-called shadow problem for v(x, t) and u(t) (cf. [10]): (2.1) Downloaded 10/23/17 to 163.1.81.103. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php For ε → 0, a one-spike steady-state solution v e , u e for (2.1), with spike centered at x = 0, is given for p > 1 by (cf. [10]) , where γ ≡ q/(p − 1), and w(y) is the unique positive homoclinic solution to To study the stability of this steady-state solution we linearize (2.1) about (2.2) by introducing v = v e +e λt Φ(x/ε) and u = u e +e λt η. After a short calculation, similar to that done in [10] and [23], but now accounting for the effect of delayed reaction kinetics, we obtain the nonlocal eigenvalue problem (NLEP) for Φ(y), with Φ → 0 as |y| → ∞, given by For the simpler case where the delayed reaction kinetics in (2.1) only arises in the inhibitor u and not in the activator v, so that v T = v in (2.1), the NLEP (2.4) is replaced by 2.1. Inhibitor delay: The spectrum near the origin with large delay. For the NLEP (2.4) we now determine the spectrum near λ = 0 in the limit of large delay T 1. We first write (2.5) as We let T → ∞ and look for eigenvalues near the origin with λ = O(T −1 ). As a result of the identity L 0 w = (p − 1)w p , we obtain from (2.6) that λ = 0 when χ s = (p − 1). Upon writing λ = c/T 1, so that µ = e −c , the condition χ s = (p − 1) determines c. In this way, we obtain the leading-order estimate for T 1 that where n = 0, ±1, ±2, . . . , and ξ > 0 is defined in (1.1b) in terms of the GM exponents (p, q, m, s). Therefore, for a fixed large value of the delay, there are many eigenvalues near the origin in the unstable right half-plane Re(λ) > 0. Downloaded 10/23/17 to 163.1.81.103. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php To determine a more refined approximation for the eigenvalues of (2.6) near λ = 0, we substitute the expansions into (2.6), expand χ s = (p − 1) + χ s1 /T + . . ., and collect terms of order O(T −1 ). This readily yields that Since L 0 w = 0, the corresponding homogeneous adjoint problem, has a nontrivial solution Ψ 1 = Ψ 1 . As such, the solvability condition for (2.9) is that Upon using (2.9) for χ s1 , the expression above can be solved for c 1 as To calculate the integral ratio I 1 /I p in (2.11) we use the readily derived identity L −1 0 w = w/(p − 1) + yw /2. We multiply (2.10) by L −1 0 w and integrate by parts, using the decay of w as y → ±∞, to obtain Then, upon using our expression for L −1 0 w, we obtain that Finally, upon integrating by parts, we get I 1 /I p = 1 − (p − 1)/(2m), which determines c 1 from (2.11). From (2.8), this yields a two-term expansion for λ in the limit of large delay T 1, For m = 2 and 1 < p ≤ 5, and with no delay T = 0, we recall from Theorem 2.3 of [23] that there is a unique critical Hopf bifurcation value τ = τ 0 H > 0 for which (2.5) has a complex conjugate pair of eigenvalues λ = ±iλ 0 IH with λ 0 IH > 0, with Re(λ) > 0 for τ > τ 0 H and Re(λ) < 0 for τ < τ 0 H . Since (2.12) shows that there are eigenvalues in Re(λ) > 0 for large delay T 1, this suggests that there should be a sequence of Hopf bifurcations at some critical values of the delay T for any τ < τ 0 H . This issue is explored in detail below for a special subrange of the GM exponents. Downloaded 10/23/17 to 163.1.81.103. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 2.2. Inhibitor delay effects: An explicitly solvable NLEP. We now show that under certain conditions on the exponents m and p the determination of the point spectrum of (2.5) in the right half-plane Re(λ) > 0 can be reduced to the study of a rather simple transcendental equation in λ. As such, for this subrange of exponents the effect of delayed inhibitor reaction kinetics is readily analyzed. We first recall some properties of the local eigenvalue problem L 0 φ l = νφ l on R for φ l ∈ H 1 (R). From [17], this problem admits the eigenvalues ν 0 > 0 and ν 1 = 0, where ν 0 is simple, and the corresponding eigenfunction φ l0 has one sign. Furthermore, from Proposition 5.6 of [4], when p ≥ 3 there are no discrete eigenvalues for L 0 in the interval (−1, 0). The continuous spectrum for L 0 is the segment λ < −1 with λ real.
As shown in Lemma 2.3 of [20], when p = 2m − 3 and m > 2, we have the key identity that We now use this identity to characterize the unstable point spectrum of the NLEP (2.5) for a particular parameter range.
Proof. We use Green's identity on w m−1 and Φ, together with the far-field decay for Φ and w as |y| → ∞, to get For eigenfunctions for which ∞ −∞ w m−1 Φ dy = 0, the first factor is nonvanishing, and we conclude To calculate the integral ratio in (2.15), we multiply w − w + w 2m−3 = 0 by w m−1 and we multiply the identity L 0 w m−1 = (m 2 − 2m)w m−1 by w. We subtract the resulting two expressions and integrate over −∞ < y < ∞. Upon using w → 0 as |y| → ∞ we obtain that (2m − 4) , which yields (2.14). Next, consider the eigenfunctions for which ∞ −∞ w m−1 Φ dy = 0. From (2.13), together with the facts that w m−1 is the unique and one-signed principal eigenfunction of L 0 , and that any eigenfunctions of the self-adjoint operator L 0 must be orthogonal, it follows that these other eigenfunctions must belong to the set of eigenfunctions of L 0 corresponding to the zero eigenvalue and any negative real eigenvalues of L 0 . Downloaded 10/23/17 to 163.1.81.103. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php We now study the roots of (2.14) in Re(λ) ≥ 0, after first writing (2.14) as (2.16) ( With no delay, so that T = 0, (2.16) reduces to the quadratic equation Since p = 2m−3, the condition (1.1b) on the GM exponents yields ξ > 0. From (2.17) we conclude that there is a unique Hopf bifurcation value τ = τ 0 H with λ = iλ 0 IH , for which Re(λ) < 0 if τ < τ 0 H and Re(λ) > 0 if τ > τ 0 H , given by Next, we show that if 0 < s < 1, there is a minimal positive value of the delay T for which (2.16) has a purely imaginary complex conjugate pair of roots when τ = 0. To show this, we set τ = 0 and λ = iλ IH in (2.16). Upon separating the real and imaginary parts of the resulting expression we obtain .
We square both sides of this expression and add, and after introducing ξ = m 2 q/(2β)− (s + 1), we derive that Since ξ > 0, and assuming that 0 < s < 1, we obtain that λ IH = λ f IH , where (2.21) Then, to determine the minimum critical value of the delay T , we substitute (2.21) into (2.19), to obtain We then use m 2 q/(2β) = ξ + (s + 1) in this expression, to obtain that T = T f , where Therefore, if 0 < s < 1, we conclude that there is a Hopf bifurcation value T f of T even when τ = 0. We remark that further Hopf bifurcations occur at the larger values Downloaded 10/23/17 to 163.1.81.103. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php of the delay T f n = T f + 2πn/λ IH for n ≥ 1, where n is an integer. For the exponent set (p, q, m, s) = (3, 2, 3, 0), for which β = 3 and ξ = 2, we get from (2.21) and (2.22) that Next, we seek a parametrization in the τ versus T plane, corresponding to the minimal value of T for which a Hopf bifurcation occurs. To do so, we write (2.16) in the form .
We now illustrate our results. For the special case (p, q, m, s) = (3, 2, 3, 0), for which ξ = 2 and β = 3, (2.25) becomes 1, we have from (2.26b) that the stability boundary satisfies τ → 0 and T → 0 as ω 0 → ∞. A plot of the stability boundary in the τ versus T plane, as computed from (2.25), is shown in the right panel of Figure 2 together with the limiting approximation given in (2.26b), which is valid for τ → 0 + . This limiting approximation is seen to provide a very accurate determination of the lower stability boundary.
In the left panel of Figure 3 we show the other Hopf bifurcations that occur for the larger values of T given by T n = T + 2πn/(3ω 0 ) for n ≥ 1 with n an integer, where T is given by (2.27). In the left panel of Figure 3, we also show the curve in the τ versus T plane where real roots of (2.16) first occur. We further remark that the large T asymptotics of (2.16), as obtained by substituting λ = c/T + c 1 /T 2 + · · · into (2.16), agrees with the general result (2.12) when p = 2m − 3. For the special case where (p, q, m, s) = (3, 2, 3, 0), (2.12) reduces to In the right panel of Figure 3, this two-term expansion for large T is favorably compared with corresponding results, computed numerically from (2.16), for τ = 0.2 and (p, q, m, s) = (3, 2, 3, 0). In summary, these results show that a time delay in the inhibitor concentration destabilizes the spike solution for the shadow problem in the sense that there exists a critical value T min of T for which the spike is unstable for all τ ≥ 0 when T > T min . Although we have only given an analysis of this for GM exponent sets where the NLEP is explicitly solvable, the numerical computations shown in Figure 1 for the standard exponent set (p, q, m, s) = (2, 1, 2, 0) reveal a similar behavior. In the next section, we show how the Hopf bifurcation boundary is readily computed numerically for arbitrary GM exponents sets for the 1-D or 2-D shadow problem. The dashed-dotted and faint dotted curves are eigenvalue paths that emerge in the right half-plane through additional Hopf bifurcation points occurring at T 1 ≈ 1.2792, and T 2 ≈ 2.4420. The black dots near the origin are the two-term asymptotic result (2.28), valid for large delay, evaluated at T = 10 for n = 0, 1, 2. N = 1, 2, assuming only a delayed reaction kinetics in the inhibitor: where |Ω| denotes the measure of Ω. A one-spike steady-state solution to (2.29) is given by (cf. [25], [10], [23]) where w(ρ) with ρ = |y| is the unique radially symmetric ground-state solution to the following BVP, in which ∆ ρ w ≡ w + (N − 1)w /ρ: (2.31) By proceeding similarly to [23], but now allowing for delayed inhibitor kinetics, the NLEP governing the stability of the steady-state solution (2.30) is It is readily shown (cf. [23]) that any unstable eigenvalue of (2.32) must be a root of g(λ) = 0, where We now seek a parametrization of the Hopf bifurcation boundary for (2.29) in the τ versus T plane. We let λ = iλ IH and set g(iλ IH ) = 0 in (2.33) to obtain By separating F(iλ IH ) into real and imaginary parts we get By taking the modulus of both sides of (2.34a), and labeling F R ≡ F R (λ IH ) and F I ≡ F I (λ IH ), we obtain that On the range of λ IH for which N s (λ IH ) > 0, we can solve for τ , and then take the imaginary part of (2.34a) to obtain the following parametric representation of the Hopf bifurcation boundary in the τ versus T plane: (2.36) From this formulation it is easy to prove that when s < 1, there is a Hopf bifurcation value T f even when τ = 0. To prove this we need only show that there is In this way, since Then, by using (2.36) for T , the Hopf bifurcation value is This establishes for our general shadow problem that whenever s < 1 there must be a Hopf bifurcation value of the delay when τ = 0 for any GM exponent set satisfying (1.1b) in either one or two dimensions. We remark that if we take p = 2m−3 for the 1-D problem where N = 1, for which the NLEP is explicitly solvable, we can readily calculate that where β = m 2 − 2m. For this special case, a simple calculation shows that the parametrization (2.36) reduces to (2.25), while (2.37b) agrees with (2.22).
We now use the parametrization (2.36) to calculate Hopf bifurcation curves in the τ versus T plane for a few GM exponent sets. To compute F R (λ IH ) and F I (λ IH ) numerically we adopt the simple numerical procedure of [23] after first computing the ground-state solution w. With this numerical procedure, the Hopf bifurcation curve for the 1-D shadow problem for the prototypical GM exponent (p, q, m, s) = (2, 1, 2, 0) is shown in Figure 1. In the left panel of Figure 4 we plot the Hopf bifurcation curve in the τ versus T parameter plane for the 2-D shadow problem for the parameter set (p, q, m, s) = (3, 2, 3, 0), which was previously considered for the 1-D case in Figure 2. With no delay, i.e., T = 0, there is a unique Hopf bifurcation value when τ ≈ 0.182, with corresponding eigenvalue λ ≈ 7.56i. In the right panel of Figure 4 we plot a similar Hopf bifurcation boundary when (p, q, m, s) = (3, 2, 3, 1) for the 2-D shadow problem. In this latter case, where s = 1, we have T → 0 as τ → 0 + , as was observed previously in Figure 3 for the 1-D shadow problem.
Next, we analyze the spectrum of the NLEP (2.32) for large delay T . Since the identity L 0 w = (p − 1)w p still holds in two dimensions, we can proceed similarly to (2.8)-(2.12) to determine the spectrum of the NLEP (2.32) for the 2-D shadow problem near the origin in the limit T 1 of large delay. We set N = 2 and readily derive that (2.39) Downloaded 10/23/17 to 163.1.81.103. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php To calculate the integral ratio I 1 /I p in (2.38) we use the identity L −1 0 w = w/(p − 1) + ρw /2 (cf. [23]), where w(ρ) satisfies (2.31). We multiply (2.39) by L −1 0 w and integrate by parts, using the decay of w as ρ → ±∞, to obtain Upon integrating the last expression by parts we get I 1 /I p = 1 − (p − 1)/m. By combining this relation with (2.38), we obtain the following two-term expression for the eigenvalues near the origin in the limit T 1 of large delay: Here c = ln(1 + ξ) + 2nπi with n = 0, ±1, ±2, . . .. Apart from a factor of 2, this expression is identical to that derived in (2.12) for the 1-D problem, and shows that many eigenvalues cluster near the origin in the limit of large delay.
3. Delay effects: An explicitly solvable NLEP problem on the infinite line. In this section we analyze the stability of a one-spike solution to the following nondimensionalized GM system with delayed reaction kinetics: , under the condition (1.1b) on the GM exponents (p, q, m, s). Since we are considering the infinite-line problem, without loss of generality we have specified a unit inhibitor diffusivity.
In [11] a one-spike steady-state solution to (3.1) was constructed for ε → 0. The result is summarized as follows.
Principal Result 3.1. For ε → 0, a one-spike steady-state solution to (3.1), labeled by v e (x) and u e (x), is given by where w(y) is the homoclinic satisfying (2.3), and where the Green's function Next, we linearize (3.1) around the steady-state solution by introducing v = v e + e λt φ and u = u e + e λt η, where φ 1 and η 1 in (3.1). This yields a singularly perturbed eigenvalue problem for φ(x) and η(x) given by In the limit ε → 0, (3.3) can be reduced to an NLEP. Since a similar calculation was done in [11] for the case of no delay, we only briefly highlight the analysis here.
Since v e is concentrated near x = 0, we look for a localized eigenfunction for φ(x) To derive the NLEP for Φ(y), we must determine η(0) from (3.3b).
We represent the solution to (3.5) in terms of G λ , and in so doing we determine η(0) as Upon substituting (3.7) into (3.4), we obtain the following NLEP for the case where both the activator and inhibitor have delayed reaction kinetics: where Φ → 0 as |y| → ∞. For the simpler case where the delayed reaction kinetics in (3.1) only arises in the inhibitor u and not in the activator v, so that v T = v in (3.1), the NLEP (3.8) is replaced by

9)
L 0 Φ ≡ Φ − Φ + pw p−1 Φ . Downloaded 10/23/17 to 163.1.81.103. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 3.1. Inhibitor delay effects: An explicitly solvable NLEP. In this subsection we analyze the spectrum of the NLEP (3.9) for the explicitly solvable case where p = 2m − 3 and m > 2. By using Lemma 2.1 we obtain that any discrete eigenvalue of the NLEP (3.9) must be a root of (3.10) For the case of no delay, where T = 0, it was shown in Principal Result 3.8 of [20] that there is a unique Hopf bifurcation value τ 0 H , with corresponding eigenvalue λ = iλ 0 IH , given by where ζ > β 2 is the smallest root of the quadratic 482. First, we determine whether there is a Hopf bifurcation value of T even when τ = 0. By setting τ = 0 in (3.10) we obtain the same problem considered in section 2.2 for the shadow problem. Therefore, for s < 1 the Hopf bifurcation value T f is given in (2.22) while the eigenvalue λ f IH is given in (2.21). For s ≥ 1, there is no Hopf bifurcation when τ = 0.
Next, we seek a parametrization of the Hopf bifurcation boundary in the τ versus T plane parametrized by λ = iλ IH , where λ IH ≡ βω 0 . We first calculate that Next, we isolate e iλ IH T in (3.10) and take the modulus of the resulting expression. By using (2.20), we obtain that Since α 2 = 1 + τ 2 β 2 ω 2 0 , we can solve for τ = τ (ω 0 ) to obtain Next, to determine T = T (ω 0 ) we separate e iλ IH T into real and imaginary parts after calculating that .
We now illustrate our result for the GM exponent set (p, q, m, s) = (3, 2, 3, 0). We obtain from (3.14) and (3.16) that In this way, in the left panel of Figure 5 we plot the Hopf bifurcation boundary in the τ versus T plane. Alternatively, for the exponent set (p, q, m, s) = (3, 2, 3, 1) for which s = 1, there is no longer any Hopf bifurcation value of T > 0 when τ = 0. For this case, a plot of the stability boundary in the τ versus T plane, as computed from (3.14) and (3.16) is shown in the right panel of Figure 5, together with the limiting approximation given in (3.17b), which is valid for τ → 0 + . For the exponent set (p, q, m, s) = (3, 2, 3, 0) in the left panel of Figure 6 we plot the other Hopf bifurcations that occur for the larger values of T given by T n = T + 2πn/(3ω 0 ) for n ≥ 1, with n an integer, where T is given in (3.18), together with the curve where real roots of (3.10) first occur. Additionally, we can fix a value of τ and calculate the path of the eigenvalues as T ranges from its first Hopf bifurcation point T 0 to some large value of T . Setting τ = 1 in (3.10), and by computing the roots of (3.10) numerically, in the right panel of Figure 6 we plot the path of some of the eigenvalues in the unstable half-plane Re(λ) > 0 as T increases. The eigenvalue path for the primary Hopf bifurcation hits the real axis at T ≈ 0.8239, with one eigenvalue tending to the origin while the other tending to λ = β = 3 as T increases. In the right panel of Figure 6 other eigenvalue branches are shown, and these branches tend to the origin as T increases.
Next, we follow the procedure in section 2.1 to determine the spectrum near the origin for the NLEP (3.9) in the limit of large delay T . The determination of the limiting asymptotics follows precisely that given in (2.7)-(2.12) except that we must replace τ by τ /2 since √ 1 + τ λ ≈ 1 + τ λ/2 + · · · for λ 1. In this way, in place of (2.12) we obtain that where c 0 ≡ ln(1 + ξ) + 2nπi with n = 0, 1, 2, . . .. For the exponent set (p, q, m, s) = (3, 2, 3, 0) of Figure 6, (3.21) reduces to Therefore, branches of eigenvalues tend to the origin as T increases, as seen in the right panel of Figure 6.

Inhibitor delay:
The infinite-line problem for arbitrary GM exponent sets. In this subsection we use the NLEP (3.9) to compute the Hopf bifurcation boundary in the τ versus T plane for a one-spike solution to the 1-D infinite-line problem (3.1) for an arbitrary exponent set (p, q, m, s) satisfying (1.1b). In analogy with (2.34), and with w(y) defined by (2.3) and L 0 Φ ≡ Φ − Φ + pw p−1 Φ, the discrete eigenvalues of the NLEP (3.9) must be the roots of (3.23) To determine the Hopf bifurcation boundaries we set λ = iλ IH in (3.23) and separate real and imaginary parts to get By taking the modulus of both sides of (3.24) we obtain that − 1 , Downloaded 10/23/17 to 163.1.81.103. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php provided that N i (λ IH ) > 0. Then, by taking the imaginary part of (3.24), restricted to the range of λ IH for which N i (λ IH ) > 0, we obtain the following parametric representation of the Hopf bifurcation boundary in the τ versus T plane: When s < 1, we can readily prove, as was done for the shadow problem studied in section 2.3, that there is a Hopf bifurcation value of T when τ = 0. To prove this we use F R (0) = 1/(p − 1) and F I (0) = 0, together with F R (λ I ) = O(λ −2 IH ) and F I (λ I ) = O(λ −1 IH ) as λ IH → ∞ (cf. [23]). This yields that N i (0) = (ξ + 1) 4 − 1 > 0 and N i → s 4 − 1 < 0 as λ IH → ∞. Therefore, there is a root λ IH = λ f IH to N i (λ IH ) = 0. Then, by using (3.26) for T , the Hopf bifurcation value of the delay is the same as given in (2.37b). Therefore, we conclude for the 1-D infinite-line problem that whenever s < 1 there must be a Hopf bifurcation value of the delay when τ = 0 for any GM exponent set satisfying (1.1b).
In the right panel of Figure 1 we use the parametrization (3.26) to calculate the Hopf bifurcation curve in the τ versus T plane for the infinite-line problem for the prototypical GM exponent set (p, q, m, s) = (2, 1, 2, 0). The numerical procedures to compute F R (λ IH ) and F I (λ IH ) are as described in section 2.3.
Finally, we remark that we can determine the spectrum of the NLEP (2.32) near the origin for large delay T . We readily derive, in analogy with (2.12) for the 1-D shadow problem, that Here c 0 = ln(1 + ξ) + 2nπi, with n = 0, ±1, ±2, . . .. This expression is identical to that derived in (2.12) for the 1-D shadow problem apart from the factor of 1/2 in τ .

Comparison with full numerics.
In order to compare our stability theory with results from full-scale numerical computations, we consider the GM model for the exponent set (p, q, m, s) = (3, 2, 3, 0) with delayed inhibitor kinetics on the finite domain |x| ≤ L, formulated as To analyze the linear stability of the one-spike steady state on an O(1) time scale we proceed as for the infinite-line problem, noting that the NLEP with this exponent set is explicitly solvable. We readily derive that the eigenvalues of the NLEP satisfy the transcendental equation Letting L → ∞ in (3.29) we recover (3.19) with m = 3 and q = 2. By setting λ = iω 0 in (3.29), and equating the real and imaginary parts of (3.29) to zero, we Downloaded 10/23/17 to 163.1.81.103. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php obtain a parametrization of the stability boundary in the τ versus T plane where a Hopf bifurcation occurs. This numerically computed boundary is shown in Figure 7 for L = 1, L = 2, and L = 500.
Next, we compare results from the stability diagram Figure 7 with full numerical results computed from (3.28). To numerically solve (3.28), we used a method of lines approach that approximated the spatial derivatives in (3.28) with centered differences and used the dde23 solver of MATLAB to solve the system of delay (ordinary) differential equations (DDEs). To generate a more precise initial condition than the asymptotic one-spike steady-state solution of Principal Result 3.1 we used this asymptotic spike profile as an initial guess for the MATLAB bvp4c solver. This numerically computed steady state was used as the initial functions on the interval t ∈ [−T, 0]. The computations shown in Figure 8 for the spike amplitude v(0, t) versus t used 151 spatial mesh points for both u and v, and the parameter set = 0.05, L = 2, D = 1, time delay T = 0.05, and the GM exponent set (p, q, m, s) = (3, 2, 3, 0). With this parameter set, the phase diagram Figure 7 predicts a Hopf bifurcation when τ = τ H ≈ 1.23, and that the spike is linearly stable only when τ < τ H . In the top left panel of Figure 8, where τ = 1.0 < τ H , we observe that the oscillation amplitude slowly decreases in time, whereas from the top right panel of Figure 8, where τ = 1.3, we observe uncontrolled temporal oscillations of the spike amplitude, followed by the collapse of the spike. In the right panel of Figure 8 corresponding to the value τ = 2.0, which is well beyond the Hopf threshold, we show a clear oscillatory collapse of the spike amplitude v(0, t). Although, as a result of a prohibitive computational expense in solving large systems of DDEs, we were unable to resolve the Hopf bifurcation threshold more precisely by increasing the number of spatial mesh points, the results in Figure 8 are entirely consistent with our linear stability predictions, and show numerically that a Hopf bifurcation occurs on the range 1 < τ < 1.3. It is an open problem to perform a weakly nonlinear analysis to determine whether the Hopf bifurcation induced by delay is in fact subcritical, as suggested by our numerical results in Figure 8 A localized spot pattern is one for which the steady-state solution for v concentrates at a discrete set of points x j ∈ Ω for j = 1, . . . , M , as → 0. In [27], the linear stability properties of such patterns were analyzed in the weak coupling regime where D = O(ν −1 ), where ν = −1/ ln . Here we extend this analysis to allow for the effect of delayed inhibitor kinetics. By setting D = D 0 /ν in (4.1), we can readily extend the analysis of [27] to show that, with delayed inhibitor kinetics, the linear stability of an M -spot solution is determined by the spectrum of the following NLEP: where ∆ ρ ≡ ∂ ρρ Φ + ρ −1 ∂ ρ Φ is the radially symmetric part of the Laplacian, and w(ρ) > 0 is the ground-state solution obtained by setting p = 2 in (2.31). In (4.2a), there are two possible choices for the multiplier χ, corresponding to either synchronous Downloaded 10/23/17 to 163.1.81.103. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php or asynchronous perturbations near each localized spot, denoted by χ s and χ a , respectively (cf. [27]): where |Ω| is the area of Ω. With either choice of the multiplier, the discrete eigenvalues of the NLEP are roots of The following lemma will be used to determine the behavior of the eigenvalues of (4.2) in the large delay limit T 1.
To consider the effect of delay, we set χ = χ a in (4.3), and let λ = iλ IH with λ IH > 0. Upon separating (4.3) into real and imaginary parts, we readily derive that at a Hopf bifurcation λ IH must be a root of (4.8a) (µ + 1) where F R (λ IH ) and F I (λ IH ) are obtained by setting N = 2 and p = m = 2 in (2.34b). Since F R (0) = 1, F I (0) = 0, while F R → 0 and F I → 0 as λ IH → ∞, the intermediate value theorem yields that there is a root to (4.8a) whenever 0 < µ < 1. As proved in section 2 of [23] we have F I (λ IH ) > 0 and F R (λ IH ) > 0 for all λ IH > 0. Moreover, from the right panel of Figure 9, where we plot the numerically computed M(λ IH ) versus λ IH , we observe that M(λ IH ) is monotone decreasing and satisfies 0 < M(λ IH ) < 1. Thus, there is no root to (4.8a), and hence no Hopf bifurcation, when µ > 1. In contrast, for any µ in 0 < µ < 1, there is a unique root to (4.8a), and the corresponding minimum value of the delay is A plot of (4.8b) versus µ on 0 < µ < 1, together with the shaded region where the spot pattern is linearly stable, is shown in the left panel of Figure 9. As a remark, since F I (λ IH ) ∼ λ IH /2 as λ IH → 0 and F R (0) = 1 (cf. [23]), we predict from (4.8b) that T ∼ 1/2 as µ → 1 − . This limiting value is indeed confirmed from the left panel of Figure 9. Next, we use (4.5) of Lemma 4.1 withχ 0 = 2/(1 + µ) andχ 1 = 0 to establish, for large delay T , that the eigenvalues near the origin have the limiting asymptotics We conclude that when µ > 1, where no Hopf bifurcation value of T occurs, we have Re(c 0 ) < 0 and so as T → ∞ the eigenvalues approach the origin from the stable left half-plane Re(λ) < 0. In contrast, when 0 < µ < 1, we have Re(c 0 ) > 0 and so, as T → ∞, the eigenvalues approach the origin from the unstable half-plane Re(λ) > 0. We note that the eigenvalue with n = 0 is real and tends to the origin as T → +∞. Finally, we determine the behavior of the eigenvalues on the positive real axis. For λ > 0 real, we write (4.3) as (4.10) (1 + µ) 2 e λT = F(λ) , Downloaded 10/23/17 to 163.1.81.103. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php where F(λ) is defined in (4.3). As proved in [23] we have that F(0) = 1, F (λ) > 0 on 0 < λ < σ 0 , where σ 0 ≈ 1.648 is the principal eigenvalue of L 0 . Moreover, it was shown in [23] that F (λ) > 0 on 0 < λ < σ 0 . We further note that the left-hand side of (4.10) begins at (1 + µ)/2, and is a monotone increasing convex function in λ for T > 0. As such, when 0 < µ < 1, there are two real roots to (4.10) on the interval 0 < λ < σ 0 when T is large enough. As T → ∞, one of these roots tends to the origin λ = 0. This root corresponds to the limiting asymptotics result (4.9) with n = 0 and T 1. In contrast, when µ > 1, there is a unique root to (4.10) on 0 < λ < σ 0 for any T > 0. This root tends to σ 0 as T → +∞.
4.2. The synchronous mode. Next, we study (4.2) with χ = χ s . From (4.3), the eigenvalues of the NLEP satisfy (4.11) e λT 2 (1 + µ) where F(λ) is defined in (4.3). To determine the Hopf bifurcation boundary we set λ = iλ IH with λ IH > 0, to obtain (4.12) e iλ IH T = 2 µ + 1 By setting the modulus of the right hand side to unity, and then solving for τ , we conclude that where |F| 2 = (F R (λ IH )) 2 + (F I (λ IH )) 2 . From the plot of |F| 2 versus λ IH shown in the right panel of Figure 9, the range of |F| 2 required in (4.13) implicitly specifies a range of λ IH . To determine the minimum value of the delay T for a Hopf bifurcation, we then separate (4.12) into real and imaginary parts. After a short calculation we obtain that (4.14) T = 1 λ IH sin −1 2 (µ + 1)(1 + τ 2 λ 2 IH ) For a fixed value of µ, the coupled system (4.13) and (4.14) specifies a curve in the τ versus T plane parametrized by λ IH . By varying µ, we then obtain a surface in the (τ, T, µ) space at which a Hopf bifurcation occurs. By using (2.34b) with N = 2 and p = m = 2 to compute F R and F I numerically, in Figure 10 we give two different views of the Hopf bifurcation surface in the (τ, T, µ) plane. We now explain some of the qualitative features of this surface.
A second qualitative feature that is evident from Figure 10 is that there is a Hopf bifurcation value of T when τ = 0 that is independent of µ. To show this analytically, we set τ = 0 in (4.12) to obtain that this threshold, denoted by T f , and the corresponding frequency λ f IH satisfy e iλ IH T = 2F R (λ IH ) + 2iF I (λ IH ), which is independent of µ. Therefore, λ IH satisfies F 2 R + F 2 I = 1/4, which yields λ f IH ≈ 2.55. The corresponding value of the delay T f = (λ f IH ) −1 tan −1 (F I /F R ) ≈ 0.379 confirms the horizontal line in the right panel of Figure 10 in the τ = 0 plane. This critical value T f separates the stability boundary surface into two sections. For 0 < T < T f , the Hopf bifurcation boundary in the τ versus µ plane has τ → τ ∞ ≈ 0.563 as µ → ∞, while τ → ∞ as µ → 1, and the region of stability is to the left of the curve. For T f < T < T max , the contour created is parabolically shaped, with τ → ∞ as µ → 0 or µ → 1, and the region of linear stability lies above the curve. In Figure 12 these differences are shown in the plots of τ versus µ for four separate values of the delay T . These plots are obtained numerically by using Newton's method on (4.13) and (4.14).  Next, we calculate the large delay behavior of the roots to the NLEP (4.2) for the synchronous mode. Sinceχ 0 = 2 andχ 1 = −2µτ /(µ + 1), we obtain from (4.5) of Lemma 4.1 that the eigenvalues near the origin have the limiting asymptotics Since Re(c 0 ) < 0, it follows for T → ∞, and for any µ > 0, that the eigenvalues approach the origin from the unstable right half-plane Re(λ) > 0. The eigenvalue in (4.17) with n = 0 is real and tends to the origin as T → +∞. Finally, we can readily show from (4.11) that, for any µ > 0, there are exactly two real eigenvalues on the interval 0 < λ < σ 0 , where σ 0 ≈ 1.648 is the principal eigenvalue of L 0 , whenever T is sufficiently large. For the exponent set (p, q, m, s) = (2, 1, 2, 0), all these different NLEPs correspond to special cases of a more general NLEP, formulated in terms of the delayed local operator where w(ρ) is the ground-state solution obtained by setting p = 2 in (2.31). In (5.1a), χ, which depends only on the product τ λ, and the dimension N , with N = 1, 2, can be chosen to correspond to one of our NLEPs as follows: In (5.1c), the two choices ofχ correspond to either the asynchronous and synchronous modes of instability for an M -spot solution in two dimensions. In this context, µ ≡ 2πM D 0 /|Ω|, as was defined in (4.2b).
It is readily shown that any unstable eigenvalue of (5.1) must be a root of g(λ) = 0, where By setting g(iλ IH ) = 0, where λ IH > 0, we obtain a 2 × 2 nonlinear system for the Hopf bifurcation values τ and λ IH at a particular value of the delay T ≥ 0. By using Newton's method on this system, in Figure 13 we plot the numerically computed Hopf bifurcation boundary in the τ versus T plane for both the 1-D shadow problem and the infinite-line problem. By comparing Figures 13 and 1, we conclude that the stability boundaries for the case of only inhibitor delay are qualitatively rather similar to those that occur when the delay arises in both the activator and inhibitor. It is intractable analytically to determine the Hopf bifurcation boundaries for the NLEP (5.1). However, under certain conditions described below, we can show that there must be a Hopf bifurcation value of the delay T . To do so, we first prove the following lemma that characterizes the behavior of certain eigenvalues of the NLEP (5.1a) for T 1.
To establish this result we proceed as follows. We first observe that eigenvalues for (5.1a) can never cross through the origin λ = 0 as either τ or T increases, owing to the fact that the eigenvalues depend on the product λT and λτ . Next, we recall from [25] that all eigenvalues of the NLEP (5.1a) satisfy Re(λ) < 0 when τ = 0 and T = 0 if and only ifχ 0 > 1. Fixing T = 0, and assumingχ 0 > 1, then, as we have shown for the different choices (5.1b) and (5.1c) in sections 2-4, (5.1a) has a Hopf bifurcation at some value τ = τ 0 H > 0, with Re(λ) < 0 for 0 < τ < τ 0 H . Finally, if we takeχ 0 > 1 and fix τ in 0 < τ < τ 0 H , then for sufficiently large delay T we have from (5.6) that many eigenvalues tend to the origin as T → ∞. Therefore, if the eigenvalue paths are continuous in T , then there must be an intermediate value of the delay T at which a Hopf bifurcation first occurred. The outline of a rigorous proof of this statement is sketched as follows.
Remark 5.1. We follow the arguments of [3] to establish the existence of Hopf bifurcation. We first check the Fredholm properties of the eigenvalue problem (5.1a).
Note that the map Φ → 2we −λT Φ−χ(τ λ)e −2λT w 2 ∞ 0 ρ N −1 wΦdρ ∞ 0 ρ N −1 w 2 dρ is a relatively compact operator of H 2 (R N ) into L 2 (R N ) since w is exponentially decaying. Hence L T − λ is Fredholm if and only if ∆ − (1 + λ) is Fredholm and this is true as long as λ > −1 or Im(λ) = 0. Thus, the operator is Fredholm of index zero if λ > −1 or Im(λ) = 0. Since |e −λT | ≤ 1 for Re(λ) ≥ 0, we conclude that all eigenvalues of (5.1a) lie on the left half-plane Re(λ) ≤ K for some K > 0, independent of T > 0. This implies the analyticity of the operator for −1 < Re(λ) ≤ K. The analyticity of the operator, the Fredholm property, and a classical theorem of Gokhberg and Krein [2,Theorem 3.6] imply that the eigenvalues of (5.1a) in the region −1 < Re(λ) < K are isolated. This then implies that the eigenvalue paths are continuous in the delay parameter T . We now choose the path containing the lowest branch in (5.6) and continue this branch in decreasing T .
Next, we remark on an open issue regarding the boundedness of the eigenvalues of the NLEP as T is increased.
Remark 5.2. It is easy to prove that if Re(λ) ≡ λ R ≥ 0, then λ R = O (1/T ) as T → ∞. In fact, if not, we would have |e −λT | → 0 which can be easily excluded by (5.1a). However, it is unclear if there holds Im(λ) ≡ λ I = O (1/T ) , as T → ∞. Indeed let λT = c T + i(b I + 2nπ), where 0 ≤ b I < 2π. Assuming for a subsequence of T n → +∞ that there holds c T → c 0 , b T → b 0 , e −λT → e −c0−ib0 = µ 0 , and 2nπ T → a 0 , then we obtain the following limiting NLEP, as T → +∞: The existence of eigenvalues of this problem would imply the existence of another branch of eigenvalues with λ = λ R + iλ I , λ R = O( 1 T ), λ I = a 0 + O( 1 T ). We leave this as an open question.
Finally, we list an open problem regarding a one-way transversal crossing condition at a Hopf bifurcation point. It is an open problem to establish whether eigenvalues of (5.1) can only enter, and never leave, the unstable right half-plane Re(λ) > 0 as the delay T increases. The establishment of such a one-way transversal crossing result that dRe(λ)/dT > 0 whenever λ = iλ IH is an eigenvalue of (5.1), combined with the fact that eigenvalue paths are continuous in T , would prove that once stability is lost at some T = T B it can never be regained at some T > T B . The numerical computations shown in the right panels of Figures 3 and 6 for the shadow and infinite-line problem, with delay only in the inhibitor kinetics, but for the different exponent set (3, 2, 3, 0), support such a one-way crossing result.
6. Discussion. Motivated by the computational studies (cf. [5], [14], [15], [18]) of pattern formation in RD systems with a time delay in the reaction kinetics, by modeling gene expression time delays, we have analyzed the linear stability of spike solutions to various limiting forms of the GM RD model with delayed reaction kinetics in both 1-D and 2-D domains. Our analysis has provided phase diagrams in parameter space where such solutions are linearly stable. We have focused our analysis on determining any Hopf bifurcation associated with the amplitude of the spike as the delay increases. Any such instability occurs on an O(1) time scale, and is associated with an unstable O(1) eigenvalue of an NLEP.
When the delay occurs only in the inhibitor kinetics, one of our main conclusions is that if the delay T exceeds a threshold value, the steady-state spike solution is unconditionally unstable. Comparison with full numerical results in one-dimension suggests that large-scale oscillations, indicative of a subcritical Hopf bifurcation, occur just beyond the Hopf bifurcation boundary. Such uncontrolled oscillations, representing a global breakdown of a robust stable patterning mechanism, were observed in Downloaded 10/23/17 to 163.1.81.103. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php the numerical studies of [5], [14], and [15] for the GM and some related models. A second main conclusion of our study is that our detailed and explicit analysis of spike stability for the special subrange of GM exponents for which the NLEP is explicitly solvable is actually representative of the more general case. More specifically, we showed that qualitatively identical conclusions regarding the stability of a spike hold for more general GM exponent sets, including the prototypical GM model where (p, q, m, s) = (2, 1, 2, 0) in (1.1). Finally, although we have primarily focused on the case where only the inhibitor kinetics has a time delay, in section 5 we have also shown that qualitatively similar results occur when both the activator and inhibitor reaction kinetics have a time delay. In particular, for this latter case where the study of the NLEP is analytically intractable, we have determined the Hopf bifurcation boundary numerically and have shown analytically, for various limiting forms of the GM model, that a Hopf bifurcation must occur as the delay increases.
We now briefly discuss a few possible extensions of this study. First, from a mathematical viewpoint, it would be interesting to investigate whether the Hopf bifurcation due to delayed reaction kinetics is typically subcritical, confirming the numerical observations in [14] and section 3.3, and therefore can lead to uncontrolled oscillations of the spike amplitude near the Hopf bifurcation boundary. Second, our analysis has been restricted to determining the linear stability of spike solutions on an O(1) time scale, as characterized by the spectrum of an NLEP. In addition, it would be interesting to study the effect of a time delay in the reaction kinetics on the small eigenvalues of order O( 2 ) (cf. [11]) in the linearization of the steady-state spike pattern. These small eigenvalues are associated with translational instabilities of the spike (cf. [11]). Any Hopf bifurcation as the delay increases that is associated with these small eigenvalues would lead to an oscillatory motion in the spatial location of the center of the spike. In a related direction, it is an open problem to derive and then analyze an ODE for the spatial location of the center of the spike for the finite-domain problem under delayed reaction kinetics. With no delay, such an analysis is given in [12].
From a modeling viewpoint, the analysis herein has shown that there is a rather restricted range of the delay T in the GM reaction kinetics that can lead to linearly stable steady-state spike patterns. Given that time delays in reaction kinetics are well motivated biologically as a result of time lags needed for gene expression (cf. [5], [14], [15], [16], [18]), a natural modeling question is how to incorporate time delays in the reaction kinetics but still maintain robustly stable spatial patterning. Although some such possible improved models are discussed in ( [5], [14], [15], [16], [18]), it would also be interesting to explore the effect of time delays on a new class of 2-D quorum-sensing models (cf. [19], [8]) for which spatially localized signaling compartments, undergoing nonlinear kinetics, are coupled through a 2-D bulk diffusion field. For a related class of 2-D models with small signaling compartments, it has been shown recently in [9] that time delays can lead to stable temporal oscillations.