A TIME-DELAY IN THE ACTIVATOR KINETICS ENHANCES THE STABILITY OF A SPIKE SOLUTION TO THE GIERER-MEINHARDT MODEL

. We study the spectrum of a new class of nonlocal eigenvalue problems (NLEPs) that characterize the linear stability properties 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 only the nonlinear autocatalytic activator kinetics. Our analysis of this model is motivated by the computational study of Seirin Lee et al. [Bull. Math. Bio., 72 (8), (2010)] on the eﬀect of gene expression time delays on spatial patterning for both the GM and some related RD models. For various limiting forms of the GM model, we show from a numerical study of the associated NLEP, together with an analytical scaling law analysis valid for large delay T , that a time-delay in only the activator kinetics is stabilizing in the sense that there is a wider region of parameter space where the spike solution is linearly stable than when there is no time delay. This enhanced stability behavior with a delayed activator kinetics is in marked contrast to the de-stabilizing eﬀect on spike solutions of having a time-delay in both the activator and inhibitor kinetics. Numerical results computed from the RD system with delayed activator kinetics are used to validate the theory for the 1-D case.


(Communicated by Shin-Ichiro Ei)
Abstract. We study the spectrum of a new class of nonlocal eigenvalue problems (NLEPs) that characterize the linear stability properties of localized spike solutions to the singularly perturbed two-component Gierer-Meinhardt (GM) reaction-diffusion (RD) system with a fixed time-delay T in only the nonlinear autocatalytic activator kinetics. Our analysis of this model is motivated by the computational study of Seirin Lee et al. [Bull. Math. Bio.,72(8), (2010)] on the effect of gene expression time delays on spatial patterning for both the GM and some related RD models. For various limiting forms of the GM model, we show from a numerical study of the associated NLEP, together with an analytical scaling law analysis valid for large delay T , that a time-delay in only the activator kinetics is stabilizing in the sense that there is a wider region of parameter space where the spike solution is linearly stable than when there is no time delay. This enhanced stability behavior with a delayed activator kinetics is in marked contrast to the de-stabilizing effect on spike solutions of having a time-delay in both the activator and inhibitor kinetics. Numerical results computed from the RD system with delayed activator kinetics are used to validate the theory for the 1-D case.

1.
Introduction. For activator-inhibitor two-component reaction-diffusion (RD) systems, it is a well-known result, originating from Turing [21], that a small perturbation of a spatially uniform steady-state solution can become unstable when the diffusivity ratio is large enough. This initial instability then leads to the generation of large-amplitude stable spatial patterns. Although this mechanism for the development of spatially inhomogeneous patterns is well-understood, and has been applied to a broad range of specific RD systems (cf. [7]) and modeling scenarios on various spatial scales, what is less well-understood is the effect on pattern development of any time-delays in the reaction kinetics. Although there are now general results for the linear stability of spatially uniform steady-states under the effect of a time-delay (cf. [2]), which have been applied to the Gierer-Meinhardt (GM) model with saturated and time-delayed reaction kinetics in [1], the effect of time-delays in the reaction kinetics on localized RD patterns is not nearly as well understood.
Time-delays in the reaction kinetics for modeling RD patterns on a cellular spatial scale are thought to be an important biological mechanism, as there typically exists a time-delay between the initiation of protein signal transduction and the time at which genes are ultimately produced (cf. [7]). In an effort to understand the effect of such gene expression time-delays on pattern formation, various new two-component activator-inhibitor RD models with a fixed time delay in the reaction kinetics were developed based on various hypothethical sub-cellular gene expression processes in [7], [15], [16], [17] (see also the survey [18]). In these studies, pattern formation aspects of the time-delayed RD models were examined through a Turing-type linear stability around a homogeneous steady-state and from large-scale numerical computations of the PDE system on both fixed and slowly growing spatial domains. For some of these activator-inhibitor models with time-delayed kinetics, related to the classical GM model and its variants (cf. [8]), it was shown that the delay induces temporal oscillations in the spatial patterning, and that these oscillations can become very large and uncontrolled as the time-delay increases, thereby suggesting a global breakdown of a robust stable patterning mechanism.
In [6], this de-stabilizing effect of a time-delay in the reaction kinetics was analyzed in detail for a general class of GM models with a time-delay in either the inhibitor kinetics or in both the activator and inhibitor kinetics. This study of [6], and the numerical simulations in [15], leads to the question of whether a time-delay in the reaction-kinetics can ever be a stabilizing effect on pattern formation. The main new goal of this paper is to show that a time-delay in only the autocatalytic term of the activator kinetics actually leads to an increase in the parameter range where localized spatial patterns are linearly stable. This stabilizing effect of a timedelay in the activator kinetics will be analyzed in detail for the GM model with prototypical exponents of the nonlinearities.
The GM model [8] for the prototypical exponent set in the semi-strong interaction limit is a singularly perturbed two-component RD system with a large diffusivity ratio. For the case where only the nonlinear activator kinetics has a finite time delay T , the activator and inhibitor concentrations v and u, respectively, in a bounded domain Ω ⊂ R N are solutions to the non-dimensional system with ∂ n u = ∂ n v = 0 for x ∈ ∂Ω. Here 0 < ε 1, v T ≡ v(x, t−T ), and the reactiontime constant τ and inhibitor diffusivity D are both O(1) positive parameters. In the limit ε → 0, it is well-known that the steady-state GM model has localized spike solutions whereby v concentrates at certain points in Ω. We will study the linear stability of these spike solutions for three different variants of this GM model with delayed activator kinetics: a one-spike solution for the 1-D shadow problem corresponding to the large inhibitor diffusivity limit D → ∞, a one-spike solution in 1-D on the infinite line, and an M -spot pattern with M ≥ 2 in a bounded 2-D domain in the weak-coupling regime D = O(− log ε). The study of these three variants of the GM model will show in a rather broad sense the stabilizing effect of delayed activator kinetics on the linear stability of localized spike solutions.
In the traditional scenario where there is no activator delay, results characterizing the existence of a Hopf bifurcation threshold value of τ have been derived for the 1-D shadow problem in [22] and [10], for the 1-D infinite-line problem in [5] and [23], and for the 2-D multi-spot problem in [26]. These previous results are all based on linearizing the RD system around a localized spike (1-D) or spot (2-D) steady-state solution and then studying the spectrum of a nonlocal eigenvalue problem (NLEP) that arises from the linearized stability problem. For the two 1-D problems, these previous results with undelayed kinetics show that a one spike-solution is linearly stable only when 0 < τ < τ H0 , and that a Hopf bifurcation occurs as τ crosses above some threshold τ H0 . Qualitatively similar results were found for the 2-D problem in [26], although the analysis and results in [26] were more intricate owing to the existence of two distinct modes of instability for the amplitudes of the spots, representing either in-phase (synchronous) or anti-phase (asynchronous) instabilities of the amplitudes of the spots. Overall, the qualitative mechanism for an oscillatory instability is that as τ increases, the inhibitor field can only respond relatively slowly to any local increase in the activator concentration due to the autocatalytic term. Such a slow response by the inhibitor field leads to an instability of the spike.
The goal of this paper is to show for each of these three variants of the GM model that the effect of a time-delay in the activator kinetics is stabilizing in the sense that there is a larger parameter range of τ , as compared with the corresponding undelayed problems, where the steady-state spike (1-D) or spot (2-D) patterns are linearly stable. In particular, for the 1-D shadow and infinite-line problems we will show, using a combined analytical and numerical study of the associated NLEP, that the Hopf bifurcation threshold for τ is a monotone increasing function of the delay T . A simple scaling law for this Hopf bifurcation threshold for the case of large delay T 1 is derived analytically for the two 1-D problems. Similar results showing that an activator time-delay has a stabilizing effect on 2-D multispot patterns are obtained from a combined analytical and numerical study of the two specific NLEPs that are associated with either asynchronous or synchronous instabilities of the amplitudes of the spots.
The mathematical challenge of this study is that the NLEP under the effect of activator delay is difficult to analyze owing to the effect of the time-delay in both the local part of the operator as well as in the multiple of the nonlocal term. In the 1-D case, this NLEP on −∞ < y < ∞ has the form for some χ(τ λ) that is analytic in Re(λ) ≥ 0. In (1.2) the delayed local operator is defined by L µ Φ ≡ Φ − Φ + 2wµΦ, where w = 3 2 sech 2 (y/2). We will provide a new detailed analytical study of spectra of the delayed local eigenvalue problem L µ φ = λφ with φ ∈ H 1 (R). By reducing this spectral problem to the study of hypergeometric functions, similar to that in [4], we derive transcendental equations characterizing all complex or real-valued spectra of L µ . After deriving a few key properties of the NLEP (1.2), and then showing how the nonlocal term eliminates unstable spectra of L µ , we use a simple numerical method together with an analytical scaling law to determine the boundary τ H = τ H (T ) in the τ versus T parameter space where (1.2) undergoes a Hopf bifurcation for both the 1-D shadow problem and 1-D infinite-line problem. Similar results for a related NLEP are derived in a 2-D context.
A key qualitative conclusion from our analysis is that a time-delay in only the activator kinetics has a stabilizing influence on the stability of localized spike solutions. This result is in marked contrast to the results found in [6] where a time-delay was assumed in both the activator and inhibitor kinetics, or only in the inhibitor kinetics. In the former case, the associated NLEP is (1.2) where χ(τ λ)µ is replaced by χ(τ λ)µ 2 , whereas in the latter case (1.2) holds, but with L µ Φ replaced with its undelayed counterpart LΦ ≡ Φ − Φ + 2wΦ. With either of these two seemingly slight modifications of the NLEP (1.2) induced by delayed inhibitor kinetics, it was shown in [6] that the Hopf bifurcation threshold τ H is monotone decreasing in the delay T , and that, moreover, there is a finite non-zero critical value T of T for which a one-spike steady-state solution in 1-D is linearly unstable for all τ ≥ 0. In this sense, [6] showed that a steady-state spike solution can become highly unstable when inhibitor delay effects are included.
The outline of this paper is as follows. In §2, we describe the three variants of the GM model with activator time-delay and, for each variant, give the specific NLEP whose spectrum characterizes the linear stability of spike solutions. In §3, we derive a few explicit and rigorous results for spectral properties of the delayed local operator and for a general class of NLEPs in 1-D and in 2-D that are associated with delayed activator kinetics. In §4, we use the NLEP to calculate the Hopf bifurcation threshold for both the 1-D shadow and infinite-line problems, and we derive a scaling law for these boundaries in the large delay limit. Hopf bifurcation thresholds associated with the 2-D NLEP, corresponding to the 2-D multi-spot problem, are studied in §5 for both synchronous and asynchronous instabilities of the spot amplitudes. A brief discussion in §6 concludes the paper.
2. Formulation of the NLEP problems. We first study the linear stability of a one-spike steady-state solution to two variants of the prototypical GM model [8] in one-space dimension when there is a time-delay in only the activator kinetics. Our first variant is the infinite-line problem where, without loss of generality, we take the inhibitor diffusivity as D = 1. This problem is formulated as . For ε → 0 it was shown in [23] that a one-spike steady-state solution v e , u e to (2.1), which is centered at x = 0, is given by where w(y) = 3 2 sech 2 (y/2) is the homoclinic solution on −∞ < y < ∞ to w − w + w 2 = 0 with w(0) > 0, w (0) = 0, and w → 0 as |y| → ∞.
In addition, we will study the so-called limiting shadow problem for v(x, t) and u(t) (cf. [10]) on the interval |x| ≤ 1, which corresponds to taking the large inhibitor diffusivity limit D → ∞ in (1.1). This shadow problem is 3) For ε → 0, a one-spike steady-state solution v e , u e for (2.3) is given by (cf. [10]) To study the linear stability of the steady-state solution for each of these two models we linearize either (2.1) or (2.3) about the steady-state by introducing v = v e + e λt Φ(x/ε), and u = u e + e λt η. After a short calculation, similar to that done in [6], we obtain that Φ(y) and λ are eigenpairs of the nonlocal eigenvalue problem (NLEP) on −∞ < y < ∞, (2.5a) where the delayed local operator, L µ , is defined by L µ Φ ≡ Φ − Φ + 2wµΦ. In (2.5a), the multiplier χ of the nonlocal term for either the infinite-line or shadow problems is In (2.5b) we must specify the principal branch of √ 1 + τ λ (cf. [23]). The NLEP (2.5) characterizes the linear stability of a one-spike solution on an O(1) time-scale to perturbations in the amplitude of the spike (cf. [23], [6]). It is readily shown that any unstable eigenvalue of (2.5), satisfying Re(λ) > 0, must be a root of g(λ) = 0, where with N = 2, we can readily extend the analysis of [26] to show that the linear stability of an M -spot solution with delayed activator kinetics, and with τ = O(1), is characterized by the spectrum of the NLEP on 0 < ρ < ∞, is the positive radially symmetric ground-state solution to ∆ ρ w − w + w 2 = 0 with w(0) > 0, w (0) = 0, and w → 0 as ρ → ∞. In (2.7a), χ(τ λ) can assume either of the two forms (cf. [26]) where β ≡ 2πM D 0 /|Ω| and |Ω| is the area of Ω. These two choices for χ correspond to either asynchronous (anti-phase) or synchronous (in-phase) instabilities of the amplitudes of the spots (cf. [26]). Both such modes of instability are possible when M ≥ 2. Although there are M − 1 possible anti-phase modes of instability of the spot amplitudes, from the leading-order asymptotic theory of [26] leading to (2.7b) with χ = 2/(1 + β), these modes have a common stability threshold. With either form for χ, the discrete eigenvalues of the NLEP (2.7) are roots g(λ) = 0, where (2.8) 3. Some general results for the delayed operator. For either N = 1 or N = 2, in this section we derive some properties for the general NLEP for Φ = Φ(y) ∈ H 1 (R N ), given by where χ(τ λ) is given in (2.5b) when N = 1 and in (2.7b) when N = 2. For simplicity of notation, for N = 1 in (3.1) we have Φ = Φ(y), w = w(y) = 3 2 sech 2 (y/2) and the integration in (3.1) is over the real line. For N = 2, we have Φ = Φ(ρ) and w = w(ρ), with ρ = |y|, and the integration in (3.1) is over R 2 . With this notation, it readily follows that the eigenvalues λ of (3.1) are the roots of g(λ) = 0, where To analyze (3.2), we must first consider some properties of the delayed local eigenvalue problem, defined by This is done in detail for N = 1 in §3.1 by using an analysis based on hypergeometric functions. Some partial results for the case N = 2 are given in §3.2.
3.1. Delayed local operator: One-dimensional case. When N = 1, we can derive transcendental equations for all of the eigenvalues (complex or real) of (3.3) by following the approach used in [25] for the case where µ is a fixed complex constant. For the convenience of the reader, the derivation of these equations is given in Appendix A. As shown in (A.10) of Appendix A (see also page 1071 of [25]), we obtain that any eigenvalue λ of (3.3) when N = 1 must be a root of one of the transcendental equations K l (λ) = 0, for l = 0, 1, 2, . . ., defined by We first observe that the translation mode λ = 0, Φ = w , must be an eigenpair for all T ≥ 0. This eigenvalue corresponds to setting l = 1 in (3.4). Next, we characterize any non-zero real-valued eigenvalue satisfying λ > −1 that exists for all T ≥ 0. It is easy to see that such eigenvalues can only occur when l = 0 or l = 2 in (3.4). In particular, for the case l = 0, we calculate that K 0 (0) < 0, K 0 (λ) → +∞ as λ → +∞, and K 0 (λ) is monotone increasing on λ > 0. Consequently, (3.4) with l = 0 has a unique root λ 0 = λ 0 (T ) > 0 for any T ≥ 0. We readily find from (3.4) that λ 0 (0) = 5/4, λ 0 ∼ log(2)/T for T 1, and λ 0 (T ) is monotone decreasing in T . In the left panel of Fig. 1 we plot λ 0 (T ) versus T , as computed numerically from (3.4). The only other real non-zero eigenvalue that exists for any T ≥ 0 is obtained from (3.4) with l = 2. Since K 2 (−1) < 0, K 2 (0) = 2 > 0, and K 2 (λ) is monotone increasing in λ, K 2 (λ) = 0 has a unique root λ 2 (T ) satisfying −1 < λ 2 (T ) < 0 for any T ≥ 0. We readily find that λ 2 → −3/4 as T → 0 and λ 2 ∼ log (3/5) /T for T 1. Next, we observe from (3.4) for l ≥ 3 that a discrete real eigenvalue emerges from the continuous spectrum λ ≤ −1 when the delay T exceeds a critical threshold T l edge > 0. By setting K l (−1) = 0 and solving for T , we identify that so that T l+1 edge > T l edge . Curiously, we find T 3 edge = 0, so that a discrete real eigenvalue emerges as soon as the delay is turned on. Since for each l ≥ 3 we have K l (0) > 0, K l (λ) is monotone increasing on −1 < λ < 0, and K l (−1) < 0 whenever T > T l edge , it follows that there is a unique root λ l (T ) to (3.4) in −1 < λ < 0 when T > T l edge . A simple calculation using (3.4) shows that, for each l ≥ 3, this eigenvalue has the following limiting asymptotics: This expression shows that all of these discrete eigenvalues that bifurcate from the continuous spectrum at critical values of the delay eventually accumulate on the stable side of the origin λ = 0 as T → ∞.   (3.10)). The path with imaginary eigenvalue when T = T j H is labeled by λ j . For even larger values of T , these paths all tend to the origin λ = 0, but in the half-space Re(λ) > 0. For each path, we also plot its continuation into Re(λ) < 0 for smaller delays.
In summary, with regards to real-valued eigenvalues of L µ when N = 1, the spectral problem (3.3) has exactly three real eigenvalues that exist for any T ≥ 0. They are λ = λ 0 (T ) > 0, λ = 0, and λ = λ 2 (T ), where −1 < λ 2 (T ) < 0. In addition, real eigenvalues bifurcate from the edge λ = −1 of the continuous spectrum at the critical values T l edge for l = 3, 4, . . . of the delay, as given explicitly in (3.5). For each l ≥ 3, this additional real eigenvalue remains in −1 < λ < 0 for all T ≥ T l edge , and it tends to the origin as T → ∞ with the asymptotic rate given in (3.6).
Next, we consider complex-valued roots of (3.4). We will only characterize those complex-valued branches of roots of (3.4) as T is varied that can exist in the unstable right-half plane Re(λ) ≥ 0. To do so, we first focus on determining any values of the delay T for which (3.3) has a pure imaginary eigenvalue λ = iλ I with λ I > 0. As shown in Appendix A, (3.4) can only have a pure imaginary root λ = iλ I for the case l = 0. As derived in Appendix A, such a pure imaginary eigenvalue occurs at the discrete values T n H , given by Here λ I , which is independent of n, is given explicitly by The uniqueness of this root follows from the fact that H(0) > 0, H(−π/3) < 0, and H(θ) is monotonic on −π/3 < θ < 0. ¿From simple numerical computations based on this explicit characterization, we obtain that θ 0 ≈ −0.99046 and λ I ≈ 2.1015. The first few critical values of the delay are (3.10) With these explicit values for which (3.4) has a pure imaginary eigenvalue when l = 0, we then readily use Newton's method on (3.4) with l = 0 and numerically path-follow branches of spectra of L µ for T > T n H . As shown in the right panel of Fig. 1 for n = 4, these branches lie in the unstable right-half plane Re(λ) > 0 and accumulate near λ = 0 as T → ∞. By path-following these same branches on the range 0 < T < T n H these branches of spectra are in Re(λ) < 0, as expected, and were created from a singular limit process T → 0 + and |λ| 1 with Re(λ) < 0. This singular behavior as T → 0 + is a well-known feature of quasipolymomials in the eigenvalue parameter that arise as the characteristic equation for traditional ODE delay equations.
Finally, we briefly discuss a general method to determine all the eigenvalues of L µ with N = 1 that exist near the origin λ = 0 in the limit T → ∞. This approach is based on the following result: Lemma 3.1. Consider the auxiliary spectral problem for Ψ(y) and ξ on −∞ < y < ∞ given by There is a countably infinite number of eigenvalues ξ l for l = 0, 1, 2, . . ., with ξ l < ξ l+1 for l ≥ 0, given explicitly by The first two eigenfunctions are Ψ 0 = w and Ψ 1 = w .
We now illustrate the use of this lemma by seeking all eigenvalues of L µ near λ = 0 when T 1. We let λ ∼ c/T for T 1 and from (3.3) obtain that To leading order we put Φ = Ψ + O(T −1 ) so that Ψ is an eigenfunction of (3.11) with eigenvalue ξ ≡ 2e −c . These eigenvalues are given in (3.12). If we set ξ 0 = 1, we have 1 = 2e −c , so that c = log(2) + 2nπi for n = 0, ±1, ±2, . . .. This gives λ ∼ [log(2) + 2nπi] /T for T 1. Setting n = 0, we obtain the asymptotics of the positive real eigenvalue λ 0 (T ) ∼ log(2)/T , derived earlier. In addition, the choices n = ±1, ±2, . . . correspond to the limiting behavior as T → ∞ of the paths of complex-valued eigenvalues in Reλ > 0 (see the right panel of Fig. 1). In contrast, if we use ξ l in (3.12) with l ≥ 2, we obtain e −c = ξ l /2 for l ≥ 2. This yields the large T behavior for the negative real eigenvalues λ l (T ) for l ≥ 2, as given in (3.6).

3.2.
Delayed local operator: Two-dimensional case. When N = 2, the explicit expression (3.4) no longer holds, and so only we must proceed indirectly and through numerical computations. For positive real eigenvalues we can still claim the following: Proof. For each µ ∈ ( 1 2 , 1), it is it is easy to see that there exists a unique principal eigenvalue, called λ(µ), to (3.3), with positive principal eigenfunction. In fact, it admits the variational characterization We only consider the range µ > 1/2, since if µ < 1/2, i.e. 2µ < 1, it holds from Lemma 5.1 (part 2) of [25] that where the equality holds iff φ = Cw. Hence for µ ≤ 1/2, λ(µ) ≤ 0 and if µ = 1/2, then λ(1/2) = 0 is the principal eigenvalue with eigenfunction w. For µ > 1/2, we have λ(µ) > 0, as can be verified by using the trial function φ = w in (3.13), which yields This implies λ(µ) > 0 when µ > 1/2. Next, for µ = 1 we recall that L 1 admits only one positive eigenvalue ν 0 , and that the second eigenvalue is zero. By a minmax variational characterization of eigenvalues based on (3.13), it follows that the eigenvalues of L µ are monotone decreasing in µ. Thus for µ < 1, the second eigenvalue L µ , is strictly less than the second eigenvalue of L 1 , which is 0. Therefore, the only positive eigenvalue is the first eigenvalue and hence it is principal, which must be simple by definition. By uniqueness, λ(µ) is a smooth function of µ.
To determine an approximation for large delay we use the scaling law ansatz λ 0 ∼ λ c /T for T 1, which yields that ∆ ρ φ − φ + 2 e −λc + · · · wφ = O(T −1 ). Therefore, φ = w + O(T −1 ), and upon using ∆ ρ w − w = −w 2 , we obtain that this equation holds when w 2 2e −λc − 1 = 0. This yields that λ c = log 2. In this way, we conclude that λ 0 (T ) has the following limiting asymptotics for small and large delay: λ 0 (T ) ∼ log 2/T + · · · , as T → ∞ . (3.17) In the left panel of Fig. 2 we show numerical results of a Newton iteration scheme applied to L µ Φ = λΦ to compute λ 0 (T ) for any T > 0 when N = 2. The limiting asyptotics (3.17) are found to compare favorably with these results. Next, we study complex eigenvalues of (3.3). We first observe that if Φ, λ is an eigenpair of (3.3) then so isΦ andλ. We first seek a necessary condition for which λ = iλ I with λ I > 0 is an eigenvalue of (3.3). We write Φ = Φ R + iΦ I and µ = e −iλ I T = µ R + iµ I in (3.3), and after separating into real and imaginary parts we get Upon multiplying the first equation by Φ I and the second by Φ R , we subtract the resulting expressions and use Green's identity to obtain Since λ I > 0, this shows that µ I = − sin(λ I T ) < 0. Next, we multiply the first equation in (3.18) by Φ R and the second by Φ I , and add the resulting equations.
so that µ R = cos(λ I T ) > 0 at any HB point. In summary, we have shown that at any HB point with λ I > 0 we must have cos(λ I T ) > 0 and sin(λ I T ) < 0. The next result establishes the existence of a HB point for L µ when N = 2. Proof. Let w(ρ) be the ground-state solution to ∆ ρ w − w + w 2 = 0 and define θ ≡ λ I T . We now show that there is a value θ 0 of θ for which has a solution λ I > 0. In fact, we consider the following eigenvalue problem where we vary θ ∈ (− π 3 , 0). When θ = 0, e −iθ = 1, it is known that (3.22) has a unique positive eigenvalue. On the other hand, when θ = −π/3, so that e −iθ = 2 , we claim that all eigenvalues of (3.22) must lie on the left-hand side of the complex λ plane. To see this, we multiply (3.22) by the conjugate of Φ and integrate to obtain the identity and − 2 sin(θ) (w|Φ| 2 ) = λ I |Φ| 2 . (3.24) Upon setting cos(θ) = 1/2, we observe from (3.23) and (3.14) that −λ R ≥ 0, so that λ R ≤ 0. If λ R = 0, we must have Φ = cw for some constant c. Substituting into the equation, we see that this is impossible. By the continuity argument (see [3]), as θ varies from 0 to − π 3 , the eigenvalue must cross the imaginary axis at some θ 0 ∈ (− π 3 , 0). From (3.23) and the fact that sin(θ) = 0 in this open interval, we have that λ = 0. Therefore the crossing point must be a Hopf bifurcation point.
From (3.24) we see that λ I > 0. Since θ 0 ∈ (− π 3 , 0), cos(θ 0 ) > 0 and sin(θ 0 ) < 0, the HB values of T are T n H ≡ (θ 0 + 2πn) /λ I for n = 1, 2, 3, . . .. ¿From simple numerical computations of a matrix eigenvalue problem obtained from a finite difference approximation of (3.21), we calculate that θ 0 ≈ −0.9303 and λ I ≈ 2.691. The first few critical values of the delay are then Paths in the complex plane for T > T j H are then similar to those shown in the right panel of Fig. 1 for the case N = 1.
We first compute that (3.26)

3.4.
Continuous dependence on T . We write (3.1) in the following form: We claim that branches of eigenvalues of this NLEP are continuous in T ≥ 0 on S ≡ {λ | Re(λ) > −1 , χ(τ λ) is analytic}. We will only work in the class of radially symmetric functions. We first look at Fredholm properties. Since the map Φ → µwΦ − χ(τ λ)w 2 is relatively compact as a map of H 2 to L 2 , we see that the operator L µ − λ is Fredholm of index zero. Next we note that, if λ 0 is an eigenvalue of L λ Φ = λΦ then its geometric multiplicity on L 2 r , where L r is the operator L µ restricted to the radial class, is 1 unless λ = λ 0 . To see this, note that if we have two linearly independent eigenfunctions, a possible combination will make wφ = 0, and hence we obtain that λ satisfies the local eigenvalue problem (3.3), which is impossible. The analyticity of our operator in λ on the set S, which can be seen from the definitions, the Fredholm property and Theorem 3.6 of Gokhberg and Krein (cf. [9]) implies that all eigenvalues of (3.31) are isolated. Finally, we show that the algebraic multiplicity is also one. This follows from Dancer's argument (see page 248 of [3]).
In conclusion, we have shown that for each fixed T , the eigenvalues of (3.31) are isolated with geometric and algebraic multiplicity one. Applying the classical theory of Kato [13], we conclude the eigenvalues vary continuously in T . To illustrate the use of this result, we let N = 1 and study numerically how the unstable complex conjugate eigenvalues for L µ are pushed into the stable left halfplane Re(λ) < 0 for the following NLEP for Φ ∈ H 1 (R 1 ) with a constant multiplier χ 0 > 0: As a starting point, we let χ 0 = 0 and choose the minimum value T 1 H of the delay for which L µ has a pure imaginary eigenvalue λ IH , and we then numerically track this HB point as χ 0 is increased. With this homotopy, the HB threshold value of T and frequency as a function of χ 0 are shown in the left and right panels of Fig. 4. From this figure, we observe that a HB occurs only when 0 ≤ χ 0 < 1, and that λ IH → 0 + and T H → +∞ as χ 0 → 1 − . For χ 0 > 1 the NLEP (3.32) does not undergo any HB as T is increased. Since when T = 0 we have Re(λ) ≤ 0, with equality holding iff Φ = w (cf. [24]), and that unstable eigenvalues can only enter Re(λ) > 0 through a HB, we conclude, by continuity in T , that Re(λ) ≤ 0 for any T > 0 whenever χ 0 > 1.
We remark that an extension of this result to N = 2 is used in §5.2 to characterize HB points for the NLEP (2.7) with the constant multiplier χ = 2/(1 + β), which applies to asynchronous instabilities of multi-spot patterns in 2-D. 4. The 1-D problems: NLEP computations and a scaling law. In this section, we consider the NLEP (2.5) for the 1-D case. From this spectral problem, we numerically determine the threshold conditions for a Hopf Bifurcation (HB) for both the infinite-line (2.1) and the shadow (2.3) problems. We also derive a scaling law for the HB thresholds in the limit of large delay T . We first note that when τ = 0 in (2.5), we have χ(τ λ) = 2, and so by the result in §3.4 (see Fig. 4), we have Re(λ) ≤ 0 for all T ≥ 0. Fixing T , we have from Theorem 3.5 that there are at least two positive real eigenvalues when τ is large enough. By continuity in τ , there must be a HB at some τ = τ H > 0 depending on T .
To determine the threshold conditions for such a HB, we set Re (g(iλ IH )) = 0 and Im (g(iλ IH )) = 0 in (2.6), where λ IH > 0. This yields a 2 × 2 nonlinear system for the HB values τ H and λ IH at a particular value of the delay T ≥ 0. By using Newton's method on this system, together with a numerical computation of (L µ − iλ IH ) −1 w 2 , in the left panel of Fig. 5 we plot the numerically computed HB boundary τ H versus T for the shadow problem (2.3). Our numerical results show that the spike solution is linearly stable for τ < τ H . The corresponding HB frequency λ IH is plotted versus T in the middle panel of Fig. 5. In the right panel of Fig. 5 we show that τ H /T and λ IH T both tend to finite non-zero limiting values as T → ∞. As shown in Fig. 6, we obtain qualitatively similar results for a HB for the infinite-line problem (2.1). Our key conclusion that τ H is monotone increasing in T shows that a time-delay in only the activator kinetics has a stabilizing effect on a spike solution for (2.3) and (2.1).
The numerical results shown in Fig. 5 and Fig. 6 show that the HB threshold and frequency satisfy τ H → ∞ and λ IH → 0 as T → ∞, respectively, with τ H /T and λ IH T both tending to finite non-zero limiting values as T → ∞.
We now characterize this large delay limiting behavior analytically. Our explicit results, described below, are shown by the dashed lines in the middle and right panels of Fig. 5 and Fig. 6. For T → ∞, we pose τ H ∼ τ 0 T and λ ∼ ic 0 /T for some c 0 > 0 and τ 0 > 0 to be found. With this scaling law, (2.5a) reduces to leading-order to where χ 0 ≡ χ(ic 0 τ 0 ). Since w − w + w 2 = 0, (4.1) yields that Φ ∼ w + O(T −1 ), provided that c 0 and τ 0 are roots to     By setting the modulus of the right-hand side of (4.3) to unity, and defining ξ ≡ c 0 τ 0 , we get 2ξ = 1 + ξ 2 , so that c 0 τ 0 = 1/ √ 3. Then, (4.3) yields e ic0 = e iπ/3 , which has c 0 = π/3 as its minimal root. In this way, the scaling law for a HB for the shadow problem (2.3) is The asymptotics for λ IH , given by the dashed curve in the middle panel of Fig. 5, agree well with the numerical results. In addition, the theoretically predicted horizontal asymptotes lim T →∞ τ H /T = τ 0 ≈ 0.551 and lim T →∞ λ IH T = c 0 ≈ 1.047, shown in the right panel of Fig. 5, agree well with the numerically computed results from (2.6).
A similar analysis can be done for the infinite-line problem where, from (2.5b), By setting the modulus of the right-hand side of (4.5) to unity, we get that 2|z −1| = |z|, where z ≡ √ 1 + iξ and ξ ≡ c 0 τ 0 . Upon separating z into real and imaginary parts, as z = z R + iz I , we get With c 0 τ 0 known, we take the imaginary part of (4.5) to obtain that sin(c 0 ) = 2z I /(z 2 R + z 2 I ). From the expressions for z R and z I in (4.6), we obtain c 0 and τ 0 in terms of α = 4(1 + √ 10)/9 as With these values of c 0 and τ 0 , the scaling law for a HB for the infinite-line problem (2.1) is τ H ∼ τ 0 T and λ = iλ IH with λ IH ∼ c 0 /T as T → ∞. The asymptotics λ IH ∼ c 0 /T is shown by the dashed curve in the middle panel of Fig. 6. The theoretically-predicted limiting values lim T →∞ τ H /T = τ 0 ≈ 1.99 and lim T →∞ λ IH T = c 0 ≈ 0.782, given by the horizontal asymptotes in the right panel of Fig. 6, agree well with the numerically computed results from (2.6).

4.1.
Numerical validation of the theory. In order to readily compare our theoretical results for the HB threshold with full numerical results from the delayed RD system we need to extend our theory to the case of a finite 1-D domain |x| ≤ L so as to be more readily able to compute solutions to the full PDE system. As such, we consider a one-spike solution centered at x = 0 for the finite-domain problem where v x = u x = 0 at x = ±L. By first constructing a one-spike solution and then analyzing the linear stability problem, a simple calculation, similar to that in [6]  , and L = 10 (heavy solid curve). The threshold was computed numerically from (2.6) with χ(τ λ) as given in (4.10). The one-spike solution is linear stable when τ < τ H . The threshold for L = 10 closely approximates that for the infinite-line problem, which was given in the left panel of Fig. 6. and [23], shows that the linear stability problem reduces to determining the roots of (2.6), where χ(τ λ) is now defined in terms of the principal branch of √ 1 + τ λ by (4.10) We then set λ = iλ IH and numerically compute the roots of g(iλ IH ) = 0 using a Newton iteration scheme, where g(λ) is defined in (2.6) in terms of χ as given in (4.10). The resulting HB curves τ H versus T for four values of L are shown in Fig. 7. As L increases we see there is a wider range of τ at a given delay T where a one-spike solution is linearly stable. As a partial verification of the results of Fig. 7 for the HB stability threshold, we took ε = 0.05, L = 2, and T = 2, and discretized (4.9) with 151 spatial meshpoints using a method of lines approach. We then used the dde23 solver of MATLAB to solve the system of delay (ordinary) differential equations (DDEs) for a value of τ slightly below and then slightly above the theoretically predicted HB threshold of τ H ≈ 5.573. The numerical results shown in the left and middle panels of Fig. 8 confirm our prediction of the HB threshold. In the right panel of Fig. 8 we show that the spike amplitude first oscillates and then collapses for a value of τ that is considerably above the HB threshold.
5. The 2-D NLEP problem: Computations and a scaling law. In this section we study the NLEP (2.7) characterizing the linear stability of an M -spot pattern, with M ≥ 2, for the GM model (1.1) with delayed activator kinetics in a bounded 2-D domain. For both the synchronous and asynchronous modes of instability, we will show that a time-delay in the activator kinetics leads to a wider parameter range where spot patterns are linearly stable. The numerics shows a slowly decaying (growing) oscillation when τ = 5.3 (τ = 5.6), respectively. A large oscillation leading to a collapse of the spike occurs when τ is well-above the HB threshold (right panel).
5.1. The synchronous mode. We first consider the NLEP (2.7) for the synchronous mode by determining a parameterization of the HB curve in the τ versus β plane for a fixed delay T ≥ 0. We substitute χ for the synchronous mode, as given in (2.7b), into (2.8) and set λ = iλ I , with λ I ≥ 0, to obtain that We solve this equation for β to get and, upon taking the imaginary part of this expression, we conclude that We then separate F µ (iλ I ) into real and imaginary parts as Upon substituting (5.4) into (5.3), we solve for τ = τ H (λ I ) to obtain the parameterization where F Rµ ≡ F Rµ (λ I ), F Iµ ≡ F Iµ (λ I ), and |F µ | ≡ F 2 Rµ + F 2 Iµ . Similarly, if we solve (5.1) for iτ λ I , we get .
(5.6) By setting the real part of the right-hand side of (5.6) to zero, and upon solving the resulting equation for β, we obtain the parameterization β = β H (λ I ) where  The expressions (5.5) and (5.7) parametrize the HB threshold in the τ versus β plane in terms of λ I > 0, at a fixed value of the delay T . We remark that if we replace F µ with F 1 , corresponding to setting T = 0 in (5.5) and (5.7), we recover the parameterization given in equation (4.16) of [6] (see also Fig. 4.1 of [6]) for the 2-D GM model with no activator or inhibitor delays.
To implement (5.5) and (5.7) numerically, we write F µ (iλ I ) in terms of the complex-valued ψ ≡ (L µ − iλ I ) −1 w 2 as where w(ρ) > 0 is the ground-state solution satisfying ∆ ρ w − w + w 2 = 0. We solve for the ground-state numerically and then determine ψ = ψ R + iψ I from a BVP solver applied to ψ. In this way, we can readily compute F Rµ (λ I ) and F Iµ (λ I ) from (5.8), which is needed in our expressions (5.5) and (5.7) for the HB threshold.
With this numerical approach, in the left panel of Fig. 9 we plot the HB threshold in the τ versus β plane for T = 0 and four nonzero values of the delay. The corresponding HB frequency, λ IH , is plotted versus β in the right panel of Fig. 9.
From this figure, we observe that a HB exists only on the range β > 1, and that τ H → +∞ as β → 1 + . We also observe that at a fixed β > 1, τ H increases as T increases, so that the effect of delayed activator kinetics is to increase the parameter range for the linear stability of the multi-spot pattern.
To determine the limiting behavior of τ H for T 1 at a fixed β > 1, we pose the scaling law τ H ∼ τ 0 T , λ ∼ ic 0 /T where c 0 > 0 and τ 0 > 0 are to be found. Upon substituting this into (3.1) we get 9) where χ 0 = χ(ic 0 τ 0 ) is given in (2.7b). Setting Φ = w + O(T −1 ), and using ∆ ρ w − w + w 2 = 0, we get −1 + (2 − χ 0 )e −ic0 = 0. This yields that Upon taking the modulus of (5.10), and solving for τ 0 c 0 , we get provided that f > 1, which implies that β > 1. Then, by taking the imaginary part of (5.10) we get In this way, we obtain for β ≡ 2πM D 0 /|Ω| > 1, that as T → +∞, there is a HB for the synchronous mode with the scaling law We observe that for a fixed T with T 1, we have τ H → +∞ and λ IH → 0 as β → 1 + . A plot of the asymptotic results for λ IH and τ H versus β on β > 1 when T = 10 are shown in Fig. 9 to compare very well with full numerical results computed from the parameterization (5.5) and (5.7).
Next, we calculate the HB threshold in the T versus β plane, obtained by solving (2.8), which we write as a coupled system for T = T H and λ IH , satisfing where F µ (iλ I ) = F Rµ (λ I ) + iF Iµ (λ I ). This system is solved as β is varied by using Newton iterations. As a starting point, we let β 1 and choose the minimum value T 1 H of the delay for which L µ has a pure imaginary eigenvalue iλ IH . We then numerically path-follow this HB point as β is decreased. With this numerical approach, the HB threshold value of T H and frequency λ IH as a function of β are shown in the left and right panels of Fig. 10. From this figure, we observe that a HB occurs only when β > 1, and that λ IH → 0 + and T H → +∞ as β → 1 + . For β < 1 the NLEP (2.7) for the asynchronous mode does not have a HB for any T ≥ 0. Although there is a HB value of T when β > 1, the entire region β > 1 is linearly unstable for any T ≥ 0, as the NLEP always has a positive real eigenvalue. In contrast, when β < 1, we conclude that Re(λ) ≤ 0 for any T ≥ 0. The HB frequency λ IH versus β. We observe that a HB occurs only on β > 1 with λ IH → 0 + and T H → +∞ as β → 1 + . For β > 1 the NLEP (2.7) for the asynchronous mode also has a positive real eigenvalue for any T ≥ 0. When β < 1, we predict that the mutli-spot pattern is linearly stable for any T ≥ 0.
Finally, we compare our HB threshold for the asynchronous mode with only activator delay with corresponding thresholds for the case of both activator and inhibitor delay, and with only inhibitor delay. In the former case v 2 T /u in (1.1) is replaced with v 2 T /u T , whereas when there is only inhibitor delay v 2 T /u and ε −N v 2 T in (1.1) are replaced with v 2 /u T and ε −N v 2 . When there is both activator and inhibitor delay, χµ in (2.7) is replaced with χµ 2 , and a HB of the NLEP must be a root λ I = λ IH and T = T H of the coupled system In contrast, when there is only inhibitor delay, L µ in (2.7) is replaced by its undelayed counterpart L 1 Φ ≡ ∆ ρ Φ − Φ + 2wΦ, and a HB of the NLEP must be a root of (1 + β) 2 cos(λ I T ) = F R1 (λ I ) , where F 1 (iλ I ) = F R1 (λ I ) + iF I1 (λ I ) now depends only on λ I and not T . As a result, in this latter case, we can conveniently parameterize the HB in the T versus β plane in terms of λ I as (5.16) By solving (5.14) using Newton's method, and by using the parameterization (5.16), in the left and right panels of Fig. 11 we plot the minimum value T H of the  Figure 11. Left panel: The minimum value T H of T versus β where a HB occurs when the time delay occurs for both the activator and inhibitor kinetics, as computed numerically from (5.14). The HB threshold occurs for any β > 0. Right panel: The corresponding HB threshold when the time-delay occurs only for the inhibitor, as computed numerically from the parameterization (5.16). The HB threshold only occurs on 0 < β < 1, and T H → 1/2 with λ IH → 0 + as β → 1 − . delay for which a HB occurs. In comparison with the left panel of Fig. 10 for the case of only activator delay, we observe from the left panel of Fig. 11 that when there is both activator and inhibitor delay there is a HB for the entire range β > 0. For the range β > 1, the HB value of T with both activator and inhibitor delay is smaller than that with only activator delay. Moreover, when there is only inhibitor delay, as shown in the right panel of Fig. 11 we observe that there is a HB only on the range 0 < β < 1, and on this range the HB threshold in T is smaller than with both activator and inhibitor delay. This HB branch terminates as β → 1 − , owing to the fact that λ IH → 0 + as β → 1 − . Since F 1I (λ I ) ∼ λ I /2 as λ I → 0 + and F R1 (0) = 1 (see [22]), we readily calculate from (5.16) that T H → 1/2 as β → 1 − . This confirms the limiting value shown in the right panel of Fig. 11.
In summary, we conclude that when there is only delayed activator kinetics the multi-spot pattern is linearly stable to asynchronous perturbations in the amplitudes of the spots for any delay T ≥ 0 when β ≡ 2πM D 0 /|Ω| < 1. If there is any inhibitor delay, then there is a HB stability threshold in T on the range 0 < β < 1. In this sense, delayed activator kinetics leads to better stability properties than when inhibitor delay is included.
6. Discussion. We have studied the onset of an oscillatory instability in the amplitude of a localized spike solution for various limiting forms of the GM activatorinhibitor RD model in the case where the nonlinear activator kinetics has a fixed time-delay. Such an instability arises from a Hopf bifurcation associated with a new class of nonlocal eigenvalue problem (NLEP). The motivation for the study of this problem is the previous computational studies (cf. [7], [15], [16], [18]) of pattern formation in RD systems with a fixed time-delay in the reaction kinetics, which models time lags needed for the expression of genes. In contrast to the conclusion of our recent analysis of [6], where a time-delay was assumed in both the activator and inhibitor kinetics, we showed herein that a time-delay in only the activator kinetics has a stabilizing effect, in the sense that there is a larger region in parameter space where a spike solution is linearly stable than when there is no delayed reaction kinetics. Phase diagrams exhibiting these larger parameter regions where the spike is linearly stable under the effect of delayed activator kinetics were generated from a numerical study of the NLEP, together with analytical scaling laws for the Hopf bifurcation thresholds in the limit of large delay.
We now briefly discuss a few possible additional directions that warrant further study. Firstly, from a modeling viewpoint, although we have considered only the effect of delayed activator kinetics on the GM model, the present study suggests more generally that, when there is a time-delay in only the autocatalytic term in the reaction kinetics, localized RD patterns will have similar enhanced stability properties. It would be interesting to study this conjecture for other choices of the reaction kinetics and, more importantly, to try to develop realistic biological modeling scenarios for which it is only the autocatalytic component of the nonlinearity that undergoes a fixed time-delay. ¿From a mathematical viewpoint, our analysis has only considered the linear stability of spike solutions on an O(1) time-scale, as characterized by the spectrum of an NLEP. For the 1-D finite-domain problem, it would be interesting to extend our analysis to study the effect of a time-delay in the activator kinetics on the small eigenvalue of order O(ε 2 ) (cf. [11]) in the linearization of a one-spike steadystate pattern. A Hopf bifurcation for this small eigenvalue as the activator delay increases would lead to a small-amplitude oscillatory motion in the spatial location of the spike. To study large-scale motion in the location of a spike, one would have to derive and then analyze a delay differential ordinary differential equation for the location of a spike. With no delayed reaction kinetics, such an analysis was given in [12].
Finally, we remark on a possible interesting effect of delayed activator kinetics on the linear stability of homoclinic stripe patterns in bounded 2-D domains. For the GM model, such homoclinic stripe solutions, formed from the localization of a spike on a one-dimensional curve in a 2-D domain, are known to be unconditionally unstable to breakup into localized spots unless one includes a strong saturation mechanism for the autocatalysis term (cf. [14], [19]). As the saturation parameter increases towards a critical value associated with a homoclinic bifurcation point, it has been shown that the principal eigenvalue of the local part of the linearized operator decreases to zero. This mechanism has been shown for straight stripes in [14], and more generally in [19], to eliminate the band of unstable breakup modes associated with the underlying NLEP, leading to a stabilization of the homoclinic stripe. Since the effect of increasing the time-delay in the activator kinetics also decreases the principal eigenvalue of the local part of the linearized operator to zero (see the left panel of Fig. 1), we anticipate that a homoclinic stripe solution will be linearly stable to breakup instabilities whenever the time-delay is large enough. It would be interesting to examine in detail this new conjectured mechanism to stabilize homoclinic stripes in 2-D domains.
Finally, to obtain λ I , as written in (3.8), we set θ = λ I T , we use 4 √ 1 + iλ I + 1 = 1 + 48(cos θ n − i sin θ n ), where θ n = −θ 0 + 2πn. Upon squaring both sides and taking the imaginary parts of the resulting expression we obtain (3.8) for λ I , which is independent of n. This completes the derivation of (3.7)-(3.9), which determines all values of the delay T for which L µ has purely imaginary eigevalues in one space dimension.