Dynamics of pulse solutions in Gierer-Meinhardt model with time dependent diffusivity Yana Nec(a) (a) Thompson Rivers University, 805 TRU Way, Kamloops, British Columbia, Canada, email: cranberryana@gmail.com Abstract Dispersive processes with a time dependent diffusivity appear in a plethora of physical systems. Most often a solution is attained for a predefined form of diffusion coefficient D(t). Here existence of pulse solutions with an arbitrary time dependence thereof is proved for the Gierer-Meinhardt model with three types of transport: regular diffusion, sub-diffusion and Lévy flights. Admission of a solution of the classical pulse shape, but for an unencumbered form of D(t) is a valuable property that allows to study phenomena of the ilk observed in various ostensibly unrelated applications. Closed form solutions are obtained for some pulse constellations. Transitions between periods of nearly constant diffusivities trigger respective cross-over between counterpart solutions known for a constant diffusivity, thereupon exhibiting otherwise unattainable behaviour, qualitatively reconstructing observable evolution peculiarities of tagged molecular structures, such as essential slowing down or speeding up during various stages of motion, inexplicable with a single constant diffusion coefficient. Keywords : non-linear pattern formation, reaction – diffusion, Gierer-Meinhardt model, pulse solution, transient diffusion regimes, anomalous diffusion 1. Background Processes with a time dependent diffusivity are used to model multifarious physical systems. Biological media often evince a non-linear dependence of the mean square displacement in time, especially in measurements of short duration (Wu and Berland, 2008). Although conventionally the diffusivity is deemed an essential quality of the diffusing substance, in complex media the influence of the environment is such that dispersive properties at a given moment in time depend on past evolution, endowing the diffusion coefficient with temporal dependence, often an algebraically decaying function D(t) ∼ t−α , 0 < α < 1. A process, whose purport is to deliver ligands, might naturally be punctuated by ingress into organelles, egress therefrom, passage through a highly crowded environment or one of a different effective dimension. An instance thereof is when molecules cross from the cytoplasm into a membrane, instantly changing a bulk diffusion, possibly in a domain of a fractional dimension, into lateral diffusion. Such transitions result in further alteration of the diffusivity function and have been accorded some attention, cf. Bakalis et al. (2015); Bronstein et al. (2009); Soula et al. (2013); Tolić-Nørrelykke et al. (2004) and Xiao et al. (2015). An entirely distinct type of process is observed in flows through porous media. There the environment’s intricate geometry renders the diffusion coefficient time dependent, whilst there exists no particular moment in time that can be designated as a cross-over and attributed to a physical change. Preprint submitted to Elsevier August 30, 2017 The diffusivity decays and saturates smoothly, and might be approximated by combinations of powers, the precise form depending on the choice of packing and other geometric parameters (Mitra et al., 1993). Similar observations are reported in Latour et al. (1993) and Latour et al. (1994). Between these two extremities lie systems with blended characteristics, cf. Cicerone et al. (1995); Luping and Gulikers (2007) and Zhang et al. (2015). What unifies all these processes is the memory, i.e. dependence of dispersive properties on past evolution. Following the cue of foregoing theoretical and experimental studies, the current contribution seeks to study transitions in a generic pattern forming system. Analysis of this ilk is expected to put forth qualitative insight into the transition mechanism and in particular the fundamental question whether based on the macroscopic measure of the mean square displacement it is possible to distinguish between transitions in the diffusant’s molecular properties and those of the medium. The system chosen must be sufficiently explored with constant D, exhibit flexibility of generalisation to diffusive regimes beyond regular Fickian diffusion and be amenable to solution with a wide class of functions D(t). In this regard the Gierer-Meinhardt model is an apposite choice. It is a chemistry based reaction – diffusion system, forming patterns as a result of the interplay between two components, activator and inhibitor (Gierer and Meinhardt, 1972). At the limit of an asymptotically small activator to inhibitor diffusivity ratio the system possesses fully non-linear solutions in the shape of localised structures known as spikes or pulses, and has been extensively studied, cf. Ni (1998), Wei (2000) and references therein. Notably this system evinces a unique diversity of possible patterns, with multiple spikes of virtually arbitrary loci (Iron and Ward, 2002). Of particular interest proved to be patterns entailing symmetry (Iron et al., 2001) and certain classes of asymmetry (Ward and Wei, 2002). Furthermore, subsequent studies analysed the system with fractional operators, corresponding to sub-diffusion (Nec and Ward, 2012) and Lévy flights (Nec, 2012). The existence of spike patterns with anomalous diffusion thus established the Gierer-Meinhardt system as a powerful paradigmatic pattern formation model: the mechanism responsible for the emergence of fully non-linear structures is so fundamental that the pattern persists even under significant modifications of the transport properties. The current contribution proves the existence of a spike pattern for the Gierer-Meinhardt system with a generic time dependent diffusion coefficient. The theory is used to exemplify the pattern evolution subject to a transition in the diffusivity and impact thereof when superimposed on an inherently varying diffusivity, as in sub-diffusion and Lévy flights types of transport. Simulation of the spike drift within the domain qualitatively reconstructs the observations referred to above of slowed down or sped up progress of tagged molecular clusters, once more attesting to the dramatic structural strength of the Gierer-Meinhardt system as a phenomenological model. The mainstay of the analysis regards the construction of a spike solution in the presence of an arbitrary diffusion coefficient function that retains the unity order of magnitude throughout the evolution process. This feature is very powerful since it gives absolute freedom to explore new types of diffusivity functions of interest if such be identified in the future. The pattern of spikes exists as a quasi-equilibrium solution, whose slow drift is governed by a differential-algebraic system of non-linear equations. §2 derives this system with regular diffusion, generalising a former result with constant diffusivity (Iron and Ward, 2002). Thereafter §3 and §4 further extend the analysis to sub-diffusive transport and Lévy flights respectively. 2 2. Spike pattern: regular diffusion Diffusion and reaction of two species, activator a(x, t) and inhibitor h(x, t) within a one dimensional finite domain is governed by ap , hq ∂t a = ǫ2 ∂x2 a − a + −1 < x < 1, t > 0 am , −1 < x < 1, t > 0, hs a(x, 0) = a0 (x), h(x, 0) = h0 (x), τo ∂t h = D(t) ∂x2 h − h + ǫ−1 ax (±1, t) = hx (±1, t) = 0, (1a) (1b) (1c) wherein the quadruple of reaction exponents (p, q, m, s) must satisfy p > 1, q > 0, m > 0, qm − s > 1, p−1 s ⩾ 0, (1d) and τo is the reaction time constant. The asymptotically small activator diffusivity ǫ2 , 0 < ǫ ≪ 1, is constant. The inhibitor diffusivity D may vary in time, however must remain of unity magnitude, i.e. D(t) ∼ O(1) ∀ t. The spike pattern is constructed by the method of matched asymptotic expansions, in itself not new (Iron and Ward, 2002), but invoked here for the much more general case of an arbitrary function D(t). The successful extension is also vital for a similar generalisation for the technically more involute analyses with sub-diffusion and Lévy flights. The outer solution for the activator is the quiescent state a(x, t) ≡ 0. For a pattern of n spikes situated  n−1 at loci xi i=0 the spatial variable of the inner layer is deduced from the form of the diffusion term in (1a), with the corresponding slow time scale reading τ = ǫ2 t. yi = (x − xi )/ǫ, (2) The inner solution is then given by  (0) (1) Ai (yi , τ ) = a xi + ǫyi , ǫ−2 τ ∼ Ai (yi , τ ) + ǫAi (yi , τ ) + . . .  (0) (1) Hi (yi , τ ) = h xi + ǫyi , ǫ−2 τ ∼ Hi (yi , τ ) + ǫHi (yi , τ ) + . . . . (3a) (3b) Substituting (3) into (1a), collecting O(1) terms and setting the boundary condition to match the outer solution, (0) p A ∂ 2 (0) (0) Ai − Ai + i(0) q = 0, 2 ∂yi Hi Similarly for the inhibitor at O ǫ−2 D τ /ǫ2   ∂ 2 (0) H = 0, ∂yi2 i lim Ai = 0. (0) (4a) (0) = lim± h(0) , (4b) |yi |−→∞ lim Hi yi −→±∞ 3 x−→xi where h(0) is the leading order term in the outer solution for the inhibitor h(x, t) ∼ h(0) (x, t) + ǫh(1) (x, t) + . . . . (4c) (0) Since D ∼ O(1) ∀ τ and Hi must be bounded as |yi | −→ ∞, equation (4b) implies that (0) Hi must be solely a function of τ . Hence it is possible to solve (4a) as (0) β (0) Ai = Hi wherein the function u u(y) = p+1 2  (τ )u yi , !1/(p−1) β= q , p−1 (5a) p−1 y 2 cosh−2/(p−1) ! (5b) is the homoclinic solution to u′′ − u + up = 0, lim u(y) = 0 (5c) |y|−→∞ and determines the spike’s shape. The spike drifts slowly, to wit the locus xi depends on τ . Correction equations read for the activator at O(ǫ) and inhibitor at O ǫ−1 respectively: ! ! (1) β dx H ∂2 du i (1) (0) − 1 + pup−1 Ai = Hi (6a) qup i(0) − ∂yi2 dyi dτ Hi D τ /ǫ2 (1) Thence ∂y2i Hi a solution  ∂ 2 (1) (0) βm−s . Hi = −um Hi 2 ∂yi (6b) must be an even function. Invoking Fredholm alternative, for (6a) to possess Z Z dxi q (1) ′2 Hi u′ up dy . (7) u dy = (0) dτ R Hi R Performing integration by parts twice for the right hand side integral, using the even property (1) of u and ∂y2i Hi , and rearranging slightly yields dxi q 1 =− dτ 2(p + 1) Hi(0) (1) (1) lim yi −→−∞ ∂Hi ∂yi + lim yi −→∞ ∂Hi ∂yi !Z up+1 dy R ,Z 2 u′ dy. (8) R The ratio of integrals a system of two linear equations Zbe computed Zby constructing Z Z in (8) can 2 up+1 dy in two unknowns, R u′ dy and R 2 u2 dy R u′ dy , by one, multiplying (5c) by R u and integrating over R, and two, multiplying (5c) by u′ and integrating twice over R. Upon solution Z Z p+1 2 p+1 , (9) u dy u′ dy = 2 p−1 R R 4 and (8) becomes q 1 dxi =− dτ p − 1 Hi(0) (1) (1) lim yi −→−∞ ∂Hi ∂yi ∂Hi ∂yi + lim yi −→∞ ! . (10) The leading order equation for the outer solution is obtained by substituting (4c) into (1b) and collecting terms of O(1): n−1 X  2 (0) (0) βm−s (0) D τ /ǫ ∂x h − h = −bm δ(x − xi ), Hi 2 bm = i=0 Z um dy, (11a) R . wherein the non-linear term ǫ−1 am hs was represented as a weighted sum of Dirac delta functions centred at the spike loci xi n−1 1 am X ∼ bi δ(x − xi ), ǫ hs i=0 1 bi = ǫ Z x+i x− i am dx = hs Z ∞ am (0) βm−s dy ∼ H i i s −∞ h Z um dy . (11b) R Equation (11a) is solved as h(0) (x, t) = n−1 X ci (τ )G(x; xi ; τ ), i=0 (0) βm−s (12a) ci = −bm Hi with G(x; xi ; τ ) being the Green’s function given by  D τ /ǫ2 G′′ − G = δ(x − xi ), G′ (±1; xi ; τ ) = 0, (12b)  an additional dependence on the slow time scale τ introduced through D τ /ǫ2 . The computation of the right hand side in (10) is completed with (1) lim yi −→−∞ −2bm    n−1 X ∂Hi ∂yi (0) βm−s Hj    j=0 j6=i (1) + lim yi −→∞ ∂Hi ∂yi = lim− ∂x h(0) + lim + ∂x h(0) = x−→xi (13) x−→xi  !   1 (0) βm−s ′ + G′ (xi ; xj ; τ ) + Hi G′ (x− . i ; xi ; τ ) + G (xi ; xi ; τ )  2   (0) Further matching of leading order inner and outer inhibitor solutions h(0) (xi , t) = Hi (τ ) in (12a) yields the differential-algebraic system for the spikes’ heights and drift of loci:     !   n−1   X βm−s βm−s 1 (0) dxi 2q bm (0) ′ ′ − ′ + G (xi ; xj ; τ ) + Hi Hj G (xi ; xi ; τ ) + G (xi ; xi ; τ ) =  dτ p − 1 Hi(0)  2    j=0  j6=i (14a) 5 (0) Hi (τ ) = −bm n−1 X (0) βm−s Hj (τ )G(xi ; xj ; τ ), j=0 i = 0, . . . n − 1. (14b) The particular case of constant diffusivity was obtained in Iron and Ward (2002). Here the dependence of G on τ is granted by D(t), a generic function throughout the solution  construction. The fact that ǫ remains present within D τ /ǫ2 in (12b) does not revoke the analysis validity since D τ /ǫ2 ∼ O(1) at all times, preserving the asymptotic orders of all expansions. When  a specific form for D(t) is prescribed, it might become more consistent to expand D τ /ǫ2 at the limit ǫ −→ 0. Nonetheless the generic form of D(t) does not interfere with an explicit solution of (12b):   √   cosh (x − 1)/ D  i  1 x+1     √ √ − x < xi cosh √    D D sinh 2/ D   G(x; xi ; τ ) = (15)   √      1 cosh (xi + 1)/ D x−1    √  √ − x > xi . cosh √    D D sinh 2/ D With (15) system (14) can be solved and integrated upon definition of a specific configuration:  n−1 number of spikes n, their initial loci xi (0) i=0 and the diffusivity D(t). At the limit ǫ −→ 0 system (14) is valid for arbitrary n and xi . In praxis, when ǫ is small, but finite, adjacent spikes must be sufficiently far apart, however since their tails decay exponentially by (5b), this requirement curbs the freedom of choice but little. Henceforth the attention is focussed on patterns with reflective symmetry, to wit ones satisfying a(x, t) = a(−x, t) and h(x, t) = h(−x, t) for all times t. Traditionally this class of spike patterns was catalogued as asymmetric, whilst symmetric patterns comprised spikes of equal heights (Ward and Wei, 2002). Patterns of reflective symmetry enable qualitative analysis of several constellations and suffice to evince non-trivial phenomena effected by transient diffusion regimes. Two forms of D(t) are considered: step discontinuous  t⩽0  0 D(t) = D1 + (D2 − D1 )θ(t − t× ), θ(t) = (16a)  1 t>0 and continuous with an inner layer D(t) = D1 + D 2 D1 − D 2 t − t× − tanh , 2 2 ε (16b) where t× is the cross-over time, and ε > 0 determines the abruptness of transition between constant inhibitor diffusivities D2 and D1 through the width of the inner layer, independently of the asymptotically small spike width ǫ. Away from the transition moment form (16a) poses a problem that is equal to a constant diffusivity case, whilst (16b) is exponentially close thereto. Form (16a) was chosen to render (14) tractable and serve as a basis for comparison due to the expected congruence of solutions with the counterparts available for a constant diffusivity, as well as for validation of numerical schemes. Form (16b) was posed to simulate 6 a more realistic smooth cross-over scenario. Its symmetry about t× is key to observations of spatial and temporal delay hereunder. In the remaining part of this section system (14) is solved in closed form with a generalisation to an arbitrary number of transitions for n = 1 and n = 2. Numerical solutions are presented for n = 3 and n = 4. 2.1. Case n = 1 Solving (14) for a single spike at an arbitrary point x0 , x0 − 1 bm x0 + 1  √  cosh √ cosh √ √ D D D sinh 2/ D (0) H0 = !−1/(βm−s−1) ,  √ . √ sinh 2x / D D 0 q dx0     . =− dτ p − 1 cosh (x − 1)/ √D cosh (x + 1)/ √D 0 0 (17a) (17b) The inference x0 = 0 is a unique stable fixed point regardless of the specific time dependence  2 D τ /ǫ follows forthwith and permits to compute the equilibrium spike height as (0) H0 = eq !1/(βm−s−1) √ 1 2 D . tanh √ bm D (18) Rearranging (17b) and integrating over the interval Iτ = [τi , τ ],  √    Z cosh 2 D √  2q dx 2x0     0′ dτ ′ = √ − D coth + τ − τi = √ p−1 dτ D sinh 2x0 D Iτ (19)  √   cosh 2 D √ 2x  √   dx , D coth √ + D sinh 2x D x0 (Iτ )  where in the latter equality D τ /ǫ2 was parameterised accordingly upon the change of integration variable. With τi = 0 without loss of generality, for the step discontinuous transition (16a) the integration can be completed analytically to yield   0 ⩽ τ ⩽ τ× F τ, 0; D1 ; 2  4qτ − (20a) =    p F τ× , 0; D1 ; 2 + F τ, τ× ; D2 ; 2 τ > τ× , Z  wherein τ× = ǫ2 t× and sinh ϑ(τ ) Dc k F τ, τi ; Dc ; k = Dc ln ln + cosh √ sinh ϑ(τi ) 2 Dc  7 cosh ϑ(τ ) − 1 cosh ϑ(τ ) + 1 ! cosh ϑ(τi ) + 1 , cosh ϑ(τi ) − 1 (20b) p with Dc being constant. For the current case of a single spike k = 2 and ϑ(τ ) = 2x0 (τ ) Dc . Expression (20a) is continuous at τ× and generalises the result for a constant diffusivity Dc −  4qτ = F τ, 0; Dc ; k p τ ⩾ 0. (21)  N Note that (20a) is readily extended for any number of step transitions N occurring at t× i i=1  N with Di i=0 the corresponding diffusivities D(t) = D0 + N   X  Di − Di−1 θ t − t× i (22a) i=1 by  i−1     P   F τ× j , τ× j−1 ; Dj−1 ; k + F τ, τ× i−1 ; Di−1 ; k    j=1 4qτ = −  p     τ× i−1 ⩽ τ ⩽ τ× i , i = 1, . . . , N     N P F τ× i , τ× i−1 ; Di−1 ; k + F τ, τ× N ; DN ; k τ > τ× N , i=1 (22b) where τ× i = ǫ t× i and τ× 0 = 0. For the continuous transition (16b) or any generic D(t) the integration in (19) must perforce be completed numerically. 2 (0) H0 H0 (0) eq,D1 (0) H0 eq,D2 x0 (τi ) x0eq x0 (τ× ) x0 D D1 D2 0 τ× τ (0) Figure 1: Single spike transition for varying diffusivity. Upper panel: spike height H0 (τ ) versus locus position x0 (τ ) for inner layer (solid) and step discontinuous (dashed) functions (equations (16b) with ε = 2 and (16a) respectively). Dotted lines show  drift towards equilibria (cross marks) for corresponding constant diffusivities D1 and D2 . Lower panel: time dependence D τ . Chemical exponents used (p, q, m, s) = (2, 1, 2, 0), with qualitatively similar results for sets (2,1,3,0), (3,2,2,0) and (4,2,2,0). Spike width ǫ = 0.1. 8 Figure 1 depicts a typical transition for a single spike. When the diffusion coefficient changes abruptly, either by a step function (16a) or continuously (16b), the spike height adjusts to the value prescribed by the corresponding constant diffusivity. A first harbinger of a transition begotten memory  is the slight delay occurring for the continuous change, (0) observable in the phase plane x0 , H0 by discerning the inequality of the area delimited by the solid and dashed lines before and after the transition point x0 (τ× ), whilst D(τ ) is symmetric about τ× . 2.2. Case n = 2 For a pair of spikes (14) prescribes two algebraic and two ordinary differential equations. To construct a pattern with reflective symmetry the spikes are set at x0 and x1 , x0 (τ ) = −x1 (τ ) (0) (0) and without loss of generality x0 < x1 . Thence (14b) gives H0 = H1 , to wit reflective symmetry entails full symmetry, and (0) H0 = whereas (14a) confirms x0 + 1 bm x  √  cosh √ 0 cosh √ √ D D D sinh 1/ D !−1/(βm−s−1) , (23a) dx1 dx0 =− with dτ dτ √  . √ sinh (2x0 + 1)/ D D dx0 q  √    . =− dτ p − 1 cosh x / D cosh (x + 1)/ √D 0 0 (23b) Here x0 = − 12 is a unique stable fixed point, and the equilibrium amplitude is (0) H0 = eq √ 2 D 1 tanh √ bm 2 D !1/(βm−s−1) . (24) For the step discontinuous transition (16a) an integration procedure similar to the quondam case of a single spike yields   F τ, 0; D1 ; 1 0 ⩽ τ ⩽ τ×  4qτ − = (25)    p F τ× , 0; D1 ; 1 + F τ, τ× ; D2 ; 1 τ > τ× , wherein the√ function F is defined in (20b) and invoked here with k = 1 and ϑ(τ ) = Dc , as before generalising (21) and readily extended to an arbitrary number 2x0 (τ ) + 1 of transitions through (22). Figure 2 depicts a typical transition for a pair of spikes. Here the delay in the response (0) Hi (x0 ) in case (16b) is well pronounced. A thorough scrutiny also reveals a faster equilibration upon transition: the solid and dashed curves reach closer to fixed points than the dotted ones for the same integration time. This further augments the asymmetry of response. The time required for the spike locus to approach the new state upon an abrupt diminution of D is surprisingly long by comparison to the rapid adjustment of height. This observation 9 x0 (0) H0 eq H0 (0) eq,D1 (0) H0 eq,D2 x0 (τ× ) x0eq x0 x0 x0 (τi ) D1 D x0 (τ× ) D2 0 x0 (τi ) τ× 0 τ× τ τ (0) Figure 2: Two spike pattern transition for varying diffusivity. Left upper panel: spike height H0 (τ ) versus locus position x0 (τ ) for inner layer (solid) and step discontinuous (dashed) functions (equations (16b) with ε = 2 and (16a) respectively). Dotted lines show drift towards equilibria (cross marks) for corresponding constant diffusivities D1 and D2 . Left lower panel: time  dependence D τ . Right panel: drift of the spike locus x0 (τ ) for inner layer (solid) and step discontinuous (dashed) functions (curves virtually indistinguishable). Dotted curves show drift with corresponding constant diffusivities. Chemical exponents used (p, q, m, s) = (2, 1, 2, 0), with qualitatively similar results for sets (2,1,3,0), (3,2,2,0) and (4,2,2,0). Spike width ǫ = 0.1. holds equally for both types of transitions. By contrast, for a single spike the four drift curves, conforming to constant diffusivities D1 and D2 and two types of transition (16), fall so closely as to be nigh indistinguishable and hence not shown. The piecewise constant extension (22) can be used as an approximation to a continuous transition. Figure 3 compares the smooth transition (16b) and piecewise constant approximation, where the span {D1 , D2 } was divided into ten equal steps, for a single spike and a pair of spikes. With 25 steps the continuous and discrete transitions are visually indistinguishable. 2.3. Case n = 3 Without loss of generality the spikes loci are at x0 < x1 < x2 . A reflectively symmetric pattern of three spikes is the first non-degenerate case, to wit possessing reflective symmetry without full symmetry, ensuing upon imposing x0 = −x2 and x1 (τ ) ≡ 0. Then (14b) for dx2 dx0 (0) (0) =− and i = 0 and i = 2 gives as expected H0 = H2 , whereas (14a) confirms dτ dτ dx1 ≡ 0. Therefore there are two non-linear algebraic equations and one ordinary differential dτ equation to be solved. Simplification of (14b) for i = 0 and i = 1 gives ! βm−s βm−s b x x + 1 m 0 0 (0) (0) (0)  √  cosh √ H0 = √ (26a) cosh √ + H1 2H0 D D 2 D sinh 1 D (0) H1 = βm−s bm 1 x0 + 1 (0) βm−s  √  2H0(0) + H1 cosh √ cosh √ √ D D 2 D sinh 1 D 10 ! . (26b) To furnish a viable numerical integration for the locus drift, some counsel is required as regards the number of roots admitted by (26) and correspondence to stability. Combining (26a) and (26b) and further simplifying yields x0 H10 cosh √ − 1 D βm−s , , H10 = −2 1 x0 + 1 H10 − cosh √ cosh √ D D (27) . (0) (0) where H10 = H1 H0 > 0 is the ratio of central to outer spike heights, by definition positive. The next concern is the, behaviour of H10 and existence of a root for (27). Note x0 1 x0 + 1 that cosh √ > 1 and cosh √ > 1 since the leftmost spike is not allowed cosh √ D D D to reach the boundary, i.e. x0 + 1 > 0 ∀ t. Thus if 0 < H10 < 1, the denominator on the right hand side.of (27) is negative, hence the numerator must be positive for a root to √  D ensues. If H10 > 1, the numerator is positive, hence exist and H10 > 1 cosh x0 , 1 x0 + 1 the denominator must be negative, giving H10 < cosh √ . For H10 = 1 cosh √ D D equation (27) simplifies to x0 + 1 2x0 + 1 = cosh √ , (28) cosh √ D D whose one solution x0 = − 32 gives the equilibrium point, and the other x0 = 0 is spurious, implying this spike merges with the central spike. In light of the above, equation (27) can be cast into , .  √  1 x0 + 1 2 H − b 10 L βm−s D, , bR = cosh √ cosh √ H10 =− , bL = 1 cosh x0 bL H10 − bR D D (29) and it becomes obvious that bL < H10 < bR . The right hand side of (29) is a convex function with a root at bL and a vertical asymptote at bR . The left hand side thereof is a convex function, as by (1d) the power exponent βm − s > 1, with a positive value at bL and positive finite value at bR . Therefore at least one root on the interval bL , bR exists. Most commonly said root is unique, however for some parameter combinations it is possible to attain three roots – one correct and the remaining pair spurious. The ratio H10 = 1 conforms to equilibrium, i.e. set of heights and loci such that no drift occurs. The drift equation (14a) with i = 0, simplified by isolation of the expressions for (0) βm−s H0 (0) βm−s and H1 from (26) and basic algebra, reads  √  ! cosh 2x0 + 1 D dx0 1 q    =− √  − H10 . dτ p − 1 √D sinh x √D cosh x0 + 1 D 0 11 (30) With this form the correspondence of x0 + 1 2x0 + 1 ≷ cosh √ cosh √ D D (31) to a pattern out of equilibrium and of (28) to the fixed point x0 = − 23 is immediate. In marked contrast to the cases n = 1 and n = 2 this fixed point is not always stable. If 0 < H10 ≶ 1 for x0 ≶ − 23 respectively at least in a close vicinity of x0 = − 23 , the point is dx0 ≷ 0. Since βm − s > 1, if 0 < H10 ≶ 1, the right hand asymptotically stable as then dτ side of (29) must satisfy the respective inequality. Thence H10 ≶ 2 + bR  . 1 + 2 bL (32) Stability ensues if this can be extended to read H10 ≶ 2 + bR  ≶ 1. 1 + 2 bL (33) Solving the rightmost inequality and rearranging, one arrives at (31), an out of equilibrium variant of (28). Using basic qualities of the hyperbolic cosine, for −1 < x0 < − 23 the top inequality in (31) holds, implying 0 < H10 < 1, i.e. central spike is smaller than the outer pair of spikes, whereas for − 23 < x0 < 0 the bottom inequality holds, conforming to H10 > 1 and central spike higher than the outer two. However it is also possible that the extension in (33) is invalid, which happens in two circumstances: 0 < H10 < 1 and the upper inequality in (32) holds, but trivially since in fact 2+b  R > 1+2 bL ; or H10 > 1 and the lower inequality in (32) holds trivially since 2 + bR < 1 + 2 bL , both entailing drift away from the fixed point. Linearising (27) by letting x0 = − 32 + δ0 , δ0 ≪ 1 and H10 = 1 + δH10 , δH10 ≪ 1 yields  . √  .√ D tanh 1 3 D 3  . √  δ0 . δH10 ∼ (34) 3 − 4(βm − s − 1) sinh2 1 3 D Although in a stable drift regime sgn δH10 = sgn δ0 , there exists a range of D such that sgn δH10 = − sgn δ0 . Completing the linearisation of (30), dδ0 ∼ −µ δ0 , dτ def µ= 6q 1 sech2 p−1 D 1 √ 3 D ! one finds that for 1 9 arcsinh2 √ βm − s − 1 !−1 1 − (βm − s − 1) sinh (35) 2  √  1/ 3 D  √  , 3 − 4(βm − s − 1) sinh2 1/ 3 D 0 is the ratio of inner to outer spike heights, by definition positive. A line of reasoning similar to the case of three spikes leads to , , x1 + 1 x1 x0 x0 + 1 . (40) bL < H10 < bR , bL = cosh √ cosh √ , bR = cosh √ cosh √ D D D D The drift equations from (14a) with i = 0 and i = 1, simplified by isolation of the expressions (0) βm−s for H0 (0) βm−s and H1 from (38) and basic algebra, read (0) βm−s−1 q bm H 0 dx0   =− dτ p − 1 D sinh 1√D (0) βm−s−1 dx1 q bm H 1   =− dτ p − 1 D sinh 1√D ! (41a) ! x1 x0 + 1 2x1 + 1 sinh √ cosh √ + sinh √ . βm−s D D D H10 (41b) x0 + 1 x1 2x0 + 1 βm−s + 2H10 sinh √ cosh √ sinh √ D D D 2 Seeking a fixed point for (41) gives βm−s −H10 eq = 2x0 + 1 sinh √ D x0 + 1 x1 2 sinh √ cosh √ D D = x1 x0 + 1 2 sinh √ cosh √ D D . 2x1 + 1 sinh √ D (42) Using the last equality and identities for hyperbolic functions one obtains x0 − x1 = 21 . Substituting this result back into (42) and combining with (39) yields , x1 + 21 x1 βm−s−1 . (43) tanh √ H10 = − tanh √ eq D D Recollecting that βm − s > 1 need not be integer, equation (43) cannot be solved in closed = 1, hence form. To circumnavigate that obstacle, it is possible to premise that H10 eq regardless of the value of βm − s equation (42) gives x1 = − 41 . This solution is then used to confirm that indeed H10 eq = 1. Thus at equilibrium the loci are at x0 = − 43 and x1 = − 14 , 14 and the heights are equal, computed from (38) as (0) (0) H0 = H1 = !1/(βm−s−1) √ 2 D 1 . tanh √ bm 4 D (44) Exploring the stability of this fixed point is slightly more involved than for patterns of fewer spikes. Taking x0 = − 43 + δ0 and x1 = − 14 + δ1 , where δ0 , δ1 ≪ 1, linearisation of (39) with H10 = 1 + δH10 , δH10 ≪ 1, yields δH10 ∼ − 1 1 δ0 + δ1 √ , tanh √ K10 4 D D 1 K10 = 2(βm − s − 1) sinh2 √ − 1 . 4 D (45) With the aid of (39) equations (41) can be recast into dx0 q bm (0) βm−s−1 = H dτ p−1 D 0 , x0 + 1 2x0 − x1 + 1 √ cosh √ D D , x0 + 1 x1 + 1 cosh √ H10 − cosh √ D D H10 − cosh q bm (0) βm−s−1 βm−s−1 dx1 = H H10 dτ p−1 D 0 2x1 − x0 √ D , x1 cosh √ D x0 1 − H10 cosh √ D , x1 cosh √ D H10 cosh (46a) . (46b) Using (45), system (46) is linearised to give 1 −3 βm−s−1   (0)  K10 q H0 d bm δ0   . ∼  δ1 dτ p − 1 D3/2 2 sinh 1 2 √D   1 +1 K10  1 +1    K10  δ0 ,   δ1 1 −3 K10  (0) (47) wherein formally H0 is the equilibrium value in (44), as taking separate variations for the (0) (0) (0) (0) two spikes’ heights Hi = Hi + δHi , δHi ≪ 1, i = {0, 1}, eq (0) (0) H1 δH10 = eq (0) −1∼ (0) (0) H0 + δH1 eq + δH0 (0) δH1 − δH0 (0) H0 . (48) eq From (47) it is apparent that the linear system has one eigenvalue that is always stable, corresponding to the matrix eigenvalue λ1 = −4,and one whose sign depends on D and 1 − 1 . The fixed point x0 = − 43 , kinetic exponents, being proportional to λ2 = 2 K10 15 x1 = − 14 , H10 = 1 is unstable if 0 < K10 < 1, corresponding to 16 arcsinh2 √ 1 βm − s − 1 !−1 0 am , −1 < x < 1, t > 0, hs a(x, 0) = a0 (x), h(x, 0) = h0 (x), τo ∂tγ h = D(t) ∂x2 h − h + ǫ−γ ax (±1, t) = hx (±1, t) = 0, Z t 1 f (t) − f (t − ζ) dγ f (t) = − dζ, γ dt Γ(−γ) 0 ζ γ+1 0<γ<1, (50a) (50b) (50c) (50d) wherein Γ(z) is the Gamma function. The fractional derivative form (50d) is often adopted in the field of dynamical systems owing to its retaining the classical property that the deriva16 tive of a constant function be the zero function, thereby permitting a generalisation of the motion near fixed points (see Oldham and Spanier (1974) for definitions that preserve other functional properties). For system (50) with a constant diffusivity Nec and Ward (2012) showed that a pattern of spikes existed, its drift near equilibrium being governed by a differential-algebraic system, whose derivation generally followed the course outlined in §2, yet was somewhat convoluted in comparison to the normal counterpart. Bearing in mind that this contribution’s purpose is to analyse the transient regime, and recalling that the construction of the pattern in §2 was successful for any function D(t) that was of order of unity, i.e. did not interfere with the orders of the terms in the asymptotic expansions, it is possible to infer forthwith an identical result for sub-diffusion. In light of the above, to extend the differential-algebraic system (7a) in Claim 2.1 of Nec and Ward (2012) to a generic function D(t) ∼ O(1) the necessary inner layer variable and slow time scale are respectively yi = x − xi , ǫγ τ = ǫγ+1 t. (51) Note that for sub-diffusion the asymptotically small activator diffusion coefficient ǫ2γ no longer equals the slow time scale order ǫγ+1 , as happens for γ = 1. This distinction is one of several consequences of the sub-diffusive transport that make system (50) a non-trivial generalisation of (1). The inner solution becomes  (0) (1) Ai (yi , τ ) = a xi + ǫγ yi , ǫ−(γ+1) τ ∼ Ai (yi , τ ) + ǫγ Ai (yi , τ ) + . . .  (0) (1) Hi (yi , τ ) = h xi + ǫγ yi , ǫ−(γ+1) τ ∼ Hi (yi , τ ) + ǫγ Hi (yi , τ ) + . . . . (52a) (52b) With these the sub-diffusive concomitant of (14) is sgn  dxi dτ  γ dxi = dτ (53a)  !   1 (0) βm−s ′ − qf (p; γ) bm (0) βm−s ′ ′ + G (xi ; xi ; τ ) + G (xi ; xi ; τ ) , G (xi ; xj ; τ ) + Hi Hj  p + 1 Hi(0)  2     j=0    n−1 X j6=i (0) Hi (τ ) = −bm n−1 X  dxi dτ (τ )G(xi ; xj ; τ ), j=0 f (p; γ) = Dγyi u ≡ sgn (0) βm−s Hj  Z 1 Γ(−γ) up+1 dyi R ,Z i = 0, . . . n − 1, u′ (yi )Dγyi u dyi , (53b) (53c) R     Z ∞ dxi dξ u(yi ) − u yi + sgn ξ . dτ ξ γ+1 0 System (50) and hence (53) inherently encode a memory feature via an effective diffusion coefficient that decays algebraically hr2 i ∼ tγ = t−α t with α = 1 − γ. A transition in the explicit coefficient D(t), entering through G(xj ; xi ; τ ), will be superimposed thereon. This 17 interaction is explored in detail hereinafter for patterns of one to four spikes. 3.1. Case n = 1 The sub-diffusive concomitant of (17) x0 + 1 x0 − 1 bm  √  cosh √ cosh √ √ D D D sinh 2/ D (0) H0 = sgn  dx0 dτ  dx0 dτ γ !−1/(βm−s−1) , (54a)  √ . √ sinh 2x / D D 0 qf (p; γ)     =− 2(p + 1) cosh (x − 1)/ √D cosh (x + 1)/ √D 0 0 (54b) cannot be integrated analytically, not even for a constant D, however the right hand side can be linearised near the equilibrium, i.e. x0 ∼ δ0 , δ0 ≪ 1, thence sgn  dδ0 dτ  dδ0 ∼ µ1/γ dτ − sgn  dδ0 dτ  δ0 !1/γ or upon integration and observation that sgn δ0 ∼      sgn δ0 (0)           δ0 (0) (γ−1)/γ δ0 (0) exp 1−γ + γ − Z τ  dδ0 dτ Z τ def qf (p; γ) , µ=  µ1/γ (τ ′ )dτ ′ µ(τ ′ )dτ ′ (55) = − sgn δ0 (0) 0 0 p+1 1 1 sech2 √ , D D ! !γ/(γ−1) 0<γ<1 (56) γ=1. For a constant D the integrals in (56) are trivial. Equation (56) allows for the inference that despite the time dependence of D the exponential and algebraic drift for normal (γ = 1) and anomalous (0 < γ < 1) transport respectively is retained for most reasonable functions D(τ ). In a simplistic case, where lim µ(τ ) = µ∞ exists (here by (16) the limit is µ∞ = D2 ), τ −→∞ it is possible to take τ > τM ≫ 1 sufficiently large, so that µ(τ ) ∼ µ∞ and write Z τM Z τ Z τM Z τ 1/γ ′ ′ 1/γ ′ ′ 1/γ ′ ′ µ (τ )dτ = µ (τ )dτ + µ (τ )dτ ∼ µ1/γ (τ ′ )dτ ′ + µ1/γ ∞ (τ − τM ) , (57) 0 τM 0 0 thence substitution into (56) reveals the time dependence of the displacement δ0 (τ ) for τ ≫ 1 as exp(−µ∞ τ ) for γ = 1 and τ γ/(γ−1) for 0 < γ < 1 (constants omitted). A more general situation, where lim µ(τ ) does not exist, will require a specific computation given τ −→∞ µ(τ ). Since µ(τ ) of a divergent magnitude or complicated singularity is unphysical, the only relevant case is that of µ(τ ) finite for all τ , but manifesting oscillatory behaviour. Depending on the particulars, a computation similar to (56) can be performed using bounds of µ(τ ), the essential functional shape of δ0 (τ ) remaining unaltered. Figure 7 shows the transition for several values of the anomaly exponent γ. By compar18 ison to γ = 1 the anomalous transition is deferred (occurs closer to equilibrium), the shift markedly non-linear. An even less trivial phenomenon is the non-monotonicity observed near the cross-over moment: juxtaposing the transition abruptness between the step discontinuous and inner layer curves, it is seen that whilst the transition is more gradual for γ = 0.25, counter-intuitively it becomes even more rapid than normal for extreme values of the anomaly index. Figure 8 compares the full numerical solution of (54) with the asymptotics (56). Evidently the asymptotic solution effectively captures the true dependence, albeit not uniformly over 0 < γ < 1: focussing the attention on the cross-over region, the only one where the error is visible, it is seen that for all values of γ the asymptotic solution predicts a slightly early transition with the least error given for normal diffusion, growing into well pronounced for significant anomaly γ ≈ 0.25, followed by improvement for an extremely deferred transition. The asymptotic solutions’ accuracy is fully correlated with the transition abruptness illustrated in figure 7. 3.2. Case n = 2 Similarly to the case of a single spike the sub-diffusive concomitant of (23) bm x0 + 1 x  √  cosh √ 0 cosh √ √ D D D sinh 1/ D (0) H0 = sgn  dx0 dτ  dx0 dτ γ !−1/(βm−s−1) ,  √ . √ sinh (2x + 1)/ D D 0 qf (p; γ)  √    =− 2(p + 1) cosh x / D cosh (x + 1)/ √D 0 0 (58a) (58b) cannot be integrated in closed form, but is amenable to a linearisation near the equilibrium x0 = − 21 + δ0 , δ0 ≪ 1, giving sgn  dδ0 dτ  dδ0 ∼ µ1/γ dτ − sgn  dδ0 dτ  δ0 !1/γ , def qf (p; γ) µ= p+1 1 1 sech2 √ , D 2 D (59) whose solution is (56) with µ as above. The delay and non-monotonic changes in transition abruptness observed for a single spike occur for a two spike pattern as well, depicted in figure 9. The asymptotic solution is indistinguishable from the full numerical solution and is not shown. As before, the approach to equilibrium is exponential for γ = 1 and algebraic for 0 < γ < 1. 3.3. Case n = 3 The equation for the ratio of central to outer spike height (27) still holds, since the system of algebraic equations (14b) is identical in form to (53b). Its solution will be different, however, due to the dependence on the spike locus, which now drifts sub-diffusively:  √  ! γ   D cosh 2x + 1 0 qf (p; γ) 1 dx0 dx0    =− sgn √  − H10 . dτ dτ 2(p + 1) √D sinh x √D cosh x + 1 D 0 0 (60) 19 Linearising (60) about x0 = − 32 and using (34) gives an equation of a form similar to that obtained for the degenerate patterns with n = 1 and n = 2: sgn  dδ0 dτ  dδ0 ∼ |µ|1/γ dτ − sgn µ sgn  dδ0 dτ  δ0 !1/γ , (61)  √  2 1 − (βm − s − 1) sinh 1/ 3 D 1 def qf (p; γ) 3  sech2 √ µ= √  , p+1 D 3 D 3 − 4(βm − s − 1) sinh2 1/ 3 D only this time including the possibility of µ < 0 that occurs for a certainrange  of D as estabdδ0 lished for regular diffusion by equation (36). Observing that sgn µ sgn = − sgn δ0 (0), dτ the solution of (61) is  !γ/(γ−1) Z τ  (γ−1)/γ 1 − γ    |µ|1/γ (τ ′ )dτ ′ 0<γ<1 + sgn µ sgn δ0 (0) δ0 (0)   γ  0 (62) δ0 ∼ !  Z  τ     µ(τ ′ )dτ ′ δ0 (0) exp − γ=1,  0 giving as before algebraic approach instead of exponential when 0 < γ < 1. Equation (60) was integrated numerically, solving (27) as required (figure 10). An interesting non-monotonicity becomes apparent: as γ diminishes, the transition occurs earlier (farther from equilibrium), however there exists a critical value of γ, below which the crossover spot is by degrees delayed. This phenomenon is explained by scrutinising equation (60) and identifying the following two opposing influences: as γ decreases, the factor f (p; γ) for a fixed value of p increases quite significantly, thereby inducing an increase in drift velocity; however, the magnitude of the right hand side remains below unity at all times, and for low values of γ the power 1/γ > 1 subdues the impact due to the growth of f (p; γ), slowing the propagation down and deferring the transition. Note the marked non-linearity of this phenomenon in γ: for γ beneath unity, but above the critical anomaly exponent (quite low and dependent on p as well as initial conditions, in the case shown approximately 0.25) the acceleration in the transition point is obvious, for γ ≈ 0.15 the cross-over for the step discontinuous change in D coincides with the normal counterpart, whereas for γ even slightly below that value the delay is manyfold more pronounced. The asymptotic solution (62) gives a good approximation only very close to the fixed point and thus is not shown. 3.4. Case n = 4 Equation (39) for the ratio of inner to outer spikes’ heights holds, anomalous transport notwithstanding. The loci of the left pair of spikes drift sub-diffusively according to sgn  dx0 dτ  γ (0) βm−s−1 dx0 qf (p; γ) bm H0   =− dτ 2(p + 1) D sinh 1√D 20 x0 + 1 x1 2x0 + 1 βm−s +2H10 sinh √ cosh √ sinh √ D D D (63a) ! sgn  dx1 dτ  γ (0) βm−s−1 dx1 qf (p; γ) bm H1   =− dτ 2(p + 1) D sinh 1√D ! x1 x0 + 1 2x1 + 1 , sinh √ cosh √ +sinh √ βm−s D D D H10 2 (63b) generalising (41). The linearised concomitant generalises (47) to  γ      1 1 dδ0 dδ0 −3 +1   sgn    K10  K10   dτ dτ   δ0   ∼µ  , γ      1  δ1 dδ1 dδ1    1 sgn +1 −3 dτ dτ K10 K10 µ= (64) 1 qf (p; γ) 1 sech2 √ . p + 1 4D 4 D With γ = 1 system (64) is linear and easily solved, both for dependence of δ0 (τ ), δ1 (τ )         1 1 δ0 1 1 −4µτ 2(1/K10 −1)µτ ∼ + (65a) δ0 (0) − δ1 (0) e δ0 (0) + δ1 (0) e δ1 −1 1 2 2  and the phase plane trajectory δ1 δ0 δ0 (τ ) + δ1 (τ ) 2K10 δ0 (τ ) − δ1 (τ ) 1−K10 = δ0 (0) + δ1 (0) 2K10 δ0 (0) − δ1 (0) 1−K10 , (65b) as long as D is constant. For a general time dependent D a numerical integration is required, albeit (65) is useful in solution comparison during the initial and post-transition stages, when D is nearly constant for the cases (16). With 0 < γ < 1 system (64) is non-linear, and in contrast to the above, allowing only for a closed form, but not explicit, solution in the phase plane with a constant D Z ( sgn v 1 + K10 + (1 − 3K10 )v 1 − 3K10 + (1 + K10 )v !1/γ  where v = δ1 δ0 . sgn v − v )−1 dv = ln δ0 + const, (66) Equation (66) was obtained by dividing the two equations in (64)  dδi and observing that the approach to equilibrium is monotonic, i.e. sgn = − sgn δi , dτ i = {0, 1}, depends on the initial point and is fixed throughout the drift, so that re-writing   γ γ−1 dδi dδi dδi dδi dδ1 sgn . The approach is mono= and division yield sgn v = sgn dτ dτ dτ dτ dδ0 dδi = 0 and vice versa for any set of parameters but K10 = 1, tonic since δi = 0 implies dτ which corresponds to the onset of instability (consult case n = 4 with regular diffusion) and not an oscillation about the equilibrium. Thus the integration in (66) is done with sgn v fixed either at unity or (−1). This is of import for the integral to be meaningful in the classical sense, similarly to the logarithm on the right hand side. Similarly to the three spike case, the linearised (fully for γ = 1 or partly for 0 < γ < 1) system (64) provides a good approximation of (63) only very near to equilibrium and for γ 21 close to normal. Figures 11 and 12 compare the deviation for γ = 1 and γ close to unity respectively. In figure 11 the linearised solution underestimates the full curve, nonetheless recreating the general form thereof correctly. This phenomenon is not transition related: observe, for instance, that for a pattern with a constant diffusivity at the initial stage the curve merges with the respective linearised solution, but not the full solution. The main purpose of this figure is to show that the linearised system does not capture the behaviour very well. This effect is further enhanced with sub-diffusion. Figure 12 depicts the evolution at the initial stage and near transition for γ = 0.95 and γ = 0.9. Whilst for γ = 0.95 the qualitative behaviour appears similar to normal, already for γ = 0.9 the partly linearised solution fails to capture the correct relative motion: the thin solid curves ascend, whereas the thin dashed curves descend. In light of the above only the full solution is used henceforward to explore the effect of sub-diffusion on the transition. Figure 13 shows the sub-diffusive counterparts of the trajectories in figure 6. As γ  diminishes, the dependence of the height ratio H10 on the loci position x0 x1 undergoes changes with respect to which spike drifts towards the equilibrium faster prior to and posttransition. For γ = 0.95 the trajectory is qualitatively similar to normal, establishing a smooth passage between γ = 1 and 0 < γ < 1 that is not trivial given the dramatic changes in the equation form with algebraic drift supplanting exponential (shown explicitly for cases  n ⩽ 3). As with regular diffusion, the outer spike drifts more slowly, so that the ratio x0 x1 increases. For γ = 0.85 at the initial stage the spikes are almost stationary whilst the height ratio first grows and subsequently diminishes. For lower γ the outer spike moves faster than the inner one, resulting in a decrease in x0 x1 during the initial stage. For γ = 0.75  the minimal value of x0 x1 is virtually equal for the step and continuous functions D(t), and for γ below   that threshold the trajectory of the continuous D(t) achieves a slightly smaller min x0 x1 . Note the trajectory for γ = 0.75 changes direction several times. The trajectories for γ = 0.65 and both transition scenarios manifest thereupon a nearly stationary period, whereas the height continues to evolve. As γ decreases further, the trajectories are qualitatively similar to that of γ = 0.65 with progressively more pronounced sluggishness of the outer spike relatively to the inner one and delayed transition for the continuous diffusion coefficient (16b) in comparison to (16a). For γ −→ 0 the integration becomes unfeasible, as f (p; γ) −→ ∞ and 1/γ −→ ∞. In summary, the main effect of a sub-diffusive transition lies in alteration of the drift velocity in all patterns explored. Sub-diffusion in its own right produces a time dependent diffusion coefficient since hr2 i ∼ tγ is equivalent to hr2 i ∼ t−α t with α = 1 − γ, i.e. effectively an algebraically decaying in time diffusivity. On this is superimposed the transition, begetting the following observations. For a single spike the cross-over is deferred, occurring closer to the equilibrium position by comparison to regular diffusion. This phenomenon is non-linear, although monotonic in γ, and accompanied by an asymmetry about the crossover point, the magnitude thereof non-monotonic in γ. Both effects are further enhanced in a constellation of two spikes. With three spikes the above shift becomes non-monotonic in γ, a counter-intuitive result. With four spikes the presence of anomaly causes reversal of the inner relative to outer spikes velocities, whilst the transition entails possibly multiple accelerations and decelerations, interspersed with quite protracted periods of virtually no motion, an occurrence otherwise unattainable. 22 4. Spike pattern: Lévy flights Lévy flights are a type of dispersion with an algebraically decaying probability density function for the particle’s step length in a continuous time random walk. A deterministic description of this transport is through the fractional Laplacian. The Gierer-Meinhardt model with the activator diffusing according to Lévy flights is ap , hq ∂t a = ǫγ Dγ|x| a − a + −1 < x < 1, t > 0, am , −1 < x < 1, t > 0, hs a(x, 0) = a0 (x), h(x, 0) = h0 (x), τo ∂t h = D(t)∂x2 h − h + ǫ−1 ax (±1, t) = hx (±1, t) = 0, a(x, t) = a(x + 2k, t), where h(x, t) = h(x + 2k, t), x ∈ R, t > 0, k ∈ Z, Z sec(πγ/2) ∞ f (x) − f (ξ) γ − dξ, D|x| f (x) = 2Γ(−γ) −∞ |x − ξ|γ+1 1<γ<2 (67a) (67b) (67c) (67d) (67e) defines the fractional Laplacian on the real line. The asymptotically small activator diffusion coefficient in (67a) ensures that at the limit ǫ −→ 0 the spike pattern exists (Nec, 2012), to wit the two species’ diffusion mechanisms retain the necessary disparity. The periodicity condition (67d) is needed for the operator (67e) to be meaningful. Similarly to the sub-diffusive case, this type of anomaly permits an immediate extension of the results for a constant diffusivity (Nec, 2012) to any D(t) ∼ O(1). Therefore with the inner solution scales (2) and expansions (3) identical to regular diffusion, the differentialalgebraic system becomes     !   n−1  dxi qf (p; γ) bm X (0) βm−s ′ 1 (0) βm−s ′ − ′ + = H G (x ; x ; τ ) + G (x ; x ; τ ) + G (x ; x ; τ ) , H i j i i i i i j  dτ p + 1 Hi(0)  2     j=0 j6=i (68a) (0) Hi (τ ) = −bm where n−1 X (0) βm−s Hj (τ )G(xi ; xj ; τ ), i = 0, . . . n − 1, j=0 f (p; γ) = Z up+1 dy R ,Z 2 u′ dy (68b) (68c) R has the same form as for the regular diffusion (9), only the function u is no longer the homoclinic solution (5b), but the least energy solution of Dγ|y| u − u + up = 0, lim u(y) = 0, |y|−→∞ (68d) introducing the dependence on γ. Thus the factor bm , defined in (11a), is now a function of γ since the spike shape u depends thereon. G(xj ; xi ; τ ) is Green’s function (15) as before. System (68) generalises (14) for Lévy flights. The single known solution to (68d) is for p = 2 23 and γ = 1, the most extreme limit of Lévy flights, u= 2 , 1 + y2 (69) whereas for any other p or 1 < γ < 2 equation (68d) must be solved numerically, precluding computation of f (p; γ) in closed form, such as (9) for the regular diffusion. Furthermore, (68d) does not accord a numerical solution as is, because the least energy solution is unstable, but can be obtained as a steady state of a system of two partial differential equations known as the Gierer-Meinhardt shadow system (consult Nec (2012) for details). Observe, however, that the factors f (p; γ) and bm (p; γ) are the sole difference between (14) and (68). Hence (68) is equivalent to (14), for instance, under a mapping of kinetic exponents (p, q, m, s) 7−→ (p̄, q̄, m̄, s̄). Bear in mind that these need not be integer, albeit classically they were so assumed. Choose p̄ and m̄ such that the factor bm remains fixed Z  Z  m m̄ bm (p, γ) = u (p,γ) dy = u (p̄,2) dy = bm̄ (p̄, 2). (70a) R R One natural choice is to keep p̄ = p, however it is possible to pick any convenient value of p̄, such as p̄ = m̄ + 1, whereby bm̄ (p̄, 2) is obtained in closed form forthwith from (5b). Further demand q f (p; γ) (p̄ + 1) qm q̄ m̄ q̄ = , s̄ = s − + . (70b) (p + 1) f (p̄; 2) p − 1 p̄ − 1 Note the resulting set (p̄, q̄, m̄, s̄) satisfies all inequalities (1d). Thus the differential-algebraic system (68), a reduction of a fractional partial differential system at the limit ǫ −→ 0, is equivalent to (14) upon a shift of chemical exponents. This is an intriguing and counter-intuitive result, given the conceptual disparity between regular Fickian diffusion and Lévy flights, where the mean square displacement is infinite. The provenance of this curious equivalence is the balance between an infinitely fast dispersion due to the operator Dγ|x| and the infinitesimally small coefficient ǫγ governing the activator’s transport in equation (67a). 5. Spectrum The linear stability analysis of patterns admitted by (1) with a constant diffusivity yields non-local eigenvalue problems. No single result encompasses the location of the spectrum for an arbitrary set of spikes due to inherent problem complexity and high dimensionality of the parameter space, but theory has been established for certain manifolds therein1 . The regular spectrum was studied in Iron et al. (2001); Ward and Wei (2002, 2003) as well as references therein, and generalised to the sub-diffusive and Lévy flights systems (50) and (67) in Nec and Ward (2012) and Nec (2012) respectively. An explicitly solvable variant is discussed in Nec and Ward (2013); Nec (2016). An inspection of these derivations reveals that the obtained non-local eigenvalue problems do not hold for an arbitrary D(t) ∼ O(1) ∀ t. 1 Consult remark 2.3 in Nec and Ward (2013). 24 By (5a) the uniformly valid quasi-equilibrium solution for the activator reads aqe = n−1 X (0) β Hj u j=0 x − x  j , (71a) G(x; xj ; τ ). (71b) ǫ and from (12a) for the inhibitor hqe = −bm n−1 X (0) βm−s Hj j=0 Define respective exponentially evolving disturbances ã and h̃ through a ∼ aqe + eλ(t) ã(x), ãx (±1, t) = 0, |ã| ≪ 1, (72a) h ∼ hqe + eλ(t) h̃(x), h̃x (±1, t) = 0, |h̃| ≪ 1. (72b) In the conventional case of constant diffusivity λ(t) = λc t with λc constant. Here this elegant simplicity proves impossible. To paraphrase, the disturbance amplitude will not be purely exponential in time, or system (1) entails no spectrum in its classical sense. Thence systems (50) and (67) will evince a similar issue. System (50) in fact possesses no spectrum even with a constant diffusivity, an aspect discussed in detail in Nec and Ward (2012), where pseudo-spectrum was defined in an apposite sense. To continue, with (72) the linear concomitant of (1) reads ′ λ (t) ã = ǫ2 ∂x2 ã − ã + p apqe ap−1 qe ã − q q+1 h̃, hqqe hqe (73a) m−1 aqe s am qe ã − h̃. s s+1 ǫ hqe ǫ hqe (73b) m τo λ (t) h̃ = D(t)∂x2 h̃ − h̃ + ′ The function ã is expected to be localised in accord with the spike constellation: n−1 x − x  X i ã ∼ ãi . ǫ i=0 (74) Changing the spatial variable in (73a) to yi = (x − xi )/ǫ and evaluating the non-linear terms (0) around the spike centre xi with the aid of (71) and hqe (xi ) = Hi gives   (0) β−1 p u h̃(xi ), i = 0, . . . , n − 1. (75a) ∂y2i ãi − 1 + λ′ (t) − pup−1 ãi = qHi Representing the non-linear compounds in (73b) as Dirac masses similarly to (11b) yields ! Z n−1   X β(m−1)−s βm−s−1 (0) (0) um−1 ãi (yi )dyi δ(x−xi ). D(t)∂x2 h̃− 1+τo λ′ (t) h̃ = h̃−mHi sbm Hi R i=0 (75b) This is a non-local drifting eigenvalue problem, where at each moment in time the growth 25 (0) rate λ′ (t) must be found, resulting in a path in the complex plain. For consistence Hi can be regarded as depending on t rather than the slow time variable τ , or equivalently D(t) and λ(t) can be expressed through τ . System (75) recovers its conventional variant by the mapping λ′ (t) 7−→ λc , thereupon rendering the eigenvalues stationary. To quantify the extent of drift it is also possible to compare the linear function λc t with the integrated λ(t), albeit to command the same ambit of detail as for the stationary problem integration is not required. Notwithstanding, for a formal completeness of the formulation an initial condition for λ(t) ought to be specified. With constant diffusivity λ(t) = λc t implies λ(0) = 0. It is convenient to adopt this condition, as then exp(λ(0)) = 1 and the disturbances defined in (72) are merely additive to the quasi-equilibria at t = 0. Any other condition would introduce a gratuitous initial scaling, but not detract from the solution validity. The drifting eigenvalue problem for the Lévy flights system (67) can be obtained from (75) by replacing the operator ∂y2i in (75a) with the fractional Laplacian (67e) and the homoclinic u by the solution of (68d). Posing the drifting pseudo-eigenvalue problem for the sub-diffusive variant (50) is not as simple and would require restrictive assumptions in addition to the limitations of the theory developed in Nec and Ward (2012). Following the analysis therein, to admit sub-diffusion expand λ(t) in (72) as λ(t) ∼ λ(0) (t) + ǫγ λ(1) (t) + · · · , 0 < γ < 1. (76) As is evident from definition (50d), a fractional derivative of an exponential is not an exponential. Therefore during the linearisation process one is required to compute the integral I = e−λ (0) (t) dγ λ(0) (t) e . dtγ (77) The drifting pseudo-eigenvalue problem will take the form (75) upon the mapping λ′ (t) 7−→ (0) (0) lim I. The computation thereof was attained for λ(0) (t) = λc t with λc constant in lemma ǫ−→0 3.1 of Nec and Ward (2012), and is not easily generalised for an arbitrary λ(0) (t). An outline is given in the appendix. The consequences of eigenvalue drift might range from minor to significant, depending on the functional form of D(t). It would be safe to conjecture that if D(t) varies reasonably slowly and remains within a previously established stability domain for a particular spike constellation at all times t, the pattern will be stable, and the opposite will hold if D(t) falls entirely in an unstable domain. When D(t) incorporates abrupt fluctuations or at some moment in time crosses from stable values to unstable or vice versa, a more complex response might be expected, including delay akin to Tzou et al. (2015) and Nec (2016), only the control parameter will be D and not τo . An exploration of the eigenvalue locations as well as ensuing dynamics is to be pursued in the future. 6. Discussion The Gierer-Meinhardt model, a paradigmatic pattern creator, was employed to provide qualitative insight into the underpinnings of a time dependent diffusion coefficient. Beginning with regular diffusion as a basic dispersion mechanism, it was shown that the spike pattern, a solution long known for the Gierer-Meinhardt model with a constant diffusivity, is in fact 26 admitted for an arbitrary diffusivity function D(t). Freedom of this ilk enables a valuable comparison in contrast to systems, where a solution is possible only for a specific form of D(t). The generalisation from a constant diffusivity to an arbitrary one was further carried out in sub-diffusive and Lévy flights concomitants of the regular system. Regular diffusion, with its linear growth of mean square displacement in time, is the apt ground for examination of the transition effect isolated from other influences. With a small number of spikes (n < 4) said effect appears unimpressive at first sight, however a closer scrutiny reveals features that could not be easily predicted, e.g. asymmetric delay in spike height adjustment for a symmetric transition and surprisingly long response time of the spike loci positions. With n = 4 the prominence of these phenomena becomes more dramatic, and in particular a markedly asymmetric evolution is obtained both for step-like and smooth cross-over. At quasi-equilibrium the variant with Lévy flights was found to be equivalent to the regular system upon re-mapping of the kinetic exponents, a non-intuitive result given that the nominal mean square displacement for Lévy flights dispersion mechanism does not exist altogether. Owing to the activator’s asymptotically small diffusion coefficient at the limit ǫ −→ 0 a proper balance is obtained nonetheless. With sub-diffusion the underlying effective diffusion coefficient decays algebraically, a type of time dependence often found in experiments, when the mean square displacement is measured over a brief period of time. A superimposed transition creates diffusivity functions requiring more than a single power for correct description, the more complicated types observed in nature. Here such combination was found to entail a non-linear as well as non-monotonic response of the evolution trajectory in the phase plane, when compared to an identical transition with regular diffusion. For instance, for the most part a waxing stronger anomaly, engendering progressively sluggish dispersion, deferred the location of transition and begot an asymmetric delay thereof for the smooth cross-over. However, with a very strong anomaly – and extremely slowly changing effective diffusion coefficient – the transition was shifted to an earlier location for a triple of spikes, whereas for a pair and a single spike the height adjustment with a smooth cross-over was nearly as simultaneous as for an abrupt cross-over. The non-linearity was best exemplified through a quick deterioration of fit between the non-linear system and partly linearised tractable versions. Another curious and arising due to the transitions phenomenon is the possibility of prolonged spells of virtual stationarity, whilst the spikes’ height nevertheless evolves substantially. During these stages the coupling between height evolution and drift is weakened. In summary, the degree of freedom accorded in an arbitrary diffusivity function revealed phenomena that are not encountered in this system with a constant diffusion coefficient. Although the studied forms of D(t) were relatively simple, the aftermath was at times counter-intuitive. The solution for a spike pattern with a general diffusivity function remains open to exploration shall other forms of D(t) be found of interest. Appendix A. Drifting pseudo-eigenvalue problem Beginning with (50d) and assuming λ(0) (t) to be sufficiently smooth to admit a power series expansion, the integral in (77) equals 1 I=− Γ(−γ) Z t( 0 1 − exp ∞ X (−1)k k=1 k! k (k) ζ λ (t) 27 !) dζ , ζ γ+1 def λ(k) (t) = dk (0) λ (t). (A.1) dtk Integrating by parts once and using power series expansion at t −→ 0+ , (" !# ∞ k X (−1) k (k) 1 1 1 − exp t λ (t) + I=− Γ(1 − γ) k! tγ k=1 ∞ X (−1)ℓ (ℓ − 1)! ℓ=1 λ(ℓ) (t) Z t ζ ℓ−γ−1 exp ∞ X (−1)k k=1 0 k! (A.2) ! ) ζ k λ(k) (t) dζ . At the limit ǫ −→ 0 the term in square brackets remains of order of unity if λ(k) are dk dk bounded, since tk k λ(0) = τ k k λ(0) . The power t−γ gives the error estimate upon replacing dt dτ   t = τ ǫγ+1 as O ǫγ(γ+1) , identically to the stationary case. Separating the first term in the def d integrand sum and changing the integration variable to ξ = λ′ (t)ζ, λ′ (t) = λ(0) (t) 6= 0, dt ∞ X (−1)ℓ 1 λ(ℓ) (t) I∼− ×  Γ(1 − γ) ℓ=1 (ℓ − 1)! λ′ (t) ℓ−γ ′ lim ǫ−→0 Z tλ (t) ξ ℓ−γ−1 exp 0 ′ " −ξ 1+ ∞ X (−1)k−1 k=2 k! λ(k) (t) ξ k−1 k λ′ (t) #!  dξ + O ǫγ(γ+1) . (A.3) The function λ (t) must be bounded. The integration path then extends infinitely far from the origin under the limiting procedure, similarly to the stationary case, only it is no longer a straight line. Due to the infinite sum in the integrand I is not as readily evaluated as for !k , dλ(0) λ(k) (t) dk λ(0) , thus O(1) for the the stationary eigenvalues. In particular k = dτ k dτ λ′ (t) entire series, unless the functional form of λ(0) is such that the higher derivatives diminish, according a possibility of an independent asymptotic expansion. This can happen only for certain premised forms of D(t). References Bakalis, E., Höfinger, S., Venturini, A., Zerbetto, F., 2015. Crossover of two power laws in the anomalous diffusion of a two lipid membrane. J. Chem. Phys. 142, 215102. Bronstein, I., Israel, Y., Kepten, E., Mai, S., Shav-Tal, Y., Barkai, E., Garini, Y., 2009. Transient anomalous diffusion of telomeres in the nucleus of mammalian cells. Phys. Rev. Lett. 103, 018102. Cicerone, M. T., Blackburn, F. R., Ediger, M. D., 1995. Anomalous diffusion of probe molecules in polystyrene: evidence for spatially heterogeneous segmental dynamics. Macromolecules 28, 8224–8232. Gierer, A., Meinhardt, H., 1972. A theory of biological pattern formation. Kybernetik 12, 30–39. 28 Iron, D., Ward, M. J., 2002. The dynamics of multi-spike solutions for the one-dimensional Gierer-Meinhardt model. SIAM J. Appl. Math. 62, 1924–1951. Iron, D., Ward, M. J., Wei, J., 2001. The stability of spike solutions to the one-dimensional Gierer-Meinhardt model. Physica D 150, 25–62. Latour, L. L., Mitra, P. P., Kleinberg, R. L., Sotak, C. H., 1993. Time-dependent diffusion coefficient of fluids in porous media as a probe of surface-to-volume ratio. J. Magnetic Resonance. Series A 101, 342–346. Latour, L. L., Svoboda, K., Mitra, P. P., Sotak, C. H., 1994. Time-dependent diffusion of water in a biological model system. Proc. Natl. Acad. Sci. USA Biophysics 91, 1229–1233. Luping, T., Gulikers, J., 2007. On the mathematics of time-dependent apparent chloride diffusion coefficient in concrete. Cement and Concrete Research 37(4), 589–595. Mitra, P. P., Sen, P. N., Schwartz, L. M., 1993. Short-time behavior of the diffusion coefficient as a geometrical probe of porous media. Phys. Rev. B 47, 8565. Nec, Y., 2012. Spike-type solutions to one dimensional Gierer-Meinhardt model with Lévy flights. Stud. Appl. Math. 129, 272–299. Nec, Y., 2016. Explicitly solvable eigenvalue problem and bifurcation delay in sub-diffusive Gierer-Meinhardt model. Euro. J. Appl. Math. DOI: 10.1017/S0956792516000012. Nec, Y., Ward, M. J., 2012. Dynamics and stability of spike-type solutions to a one dimensional Gierer-Meinhardt model with sub-diffusion. Physica D 241, 947–963. Nec, Y., Ward, M. J., 2013. An explicitly solvable nonlocal eigenvalue problem and the stability of a spike for a sub-diffusive reaction-diffusion system. Math. Model. Nat. Phenom. 8, 55–87. Ni, W., 1998. Diffusion, cross-diffusion and their spike-layer steady states. Notices of AMS 45, 9–18. Oldham, K. B., Spanier, J., 1974. Fractional Calculus. Academic Press, New York. Soula, H., Caré, B., Beslon, G., Berry, H., 2013. Anomalous versus slowed-down Brownian diffusion in the ligand-binding equilibrium. Biophys. J. 105, 2064–2073. Tolić-Nørrelykke, I. M., Munteanu, E., Thon, G., Oddershede, L., Berg-Sørensen, K., 2004. Anomalous diffusion in living yeast cells. Phys. Rev. Lett. 93, 078102. Tzou, J. C., Ward, M. J., Kolokolnikov, T., 2015. Slowly varying control parameters, delayed bifurcations, and the stability of spikes in reaction-diffusion systems. Physica D 290, 24–43. Ward, M. J., Wei, J., 2002. Asymmetric spike patterns for the one-dimensional GiererMeinhardt model: equilibria and stability. Euro. J. Appl. Math. 13, 283–320. Ward, M. J., Wei, J., 2003. Hopf bifurcation and oscillatory instabilities of spike solutions for the one-dimensional Gierer-Meinhardt model. J. Nonl. Sci. 13, 209–264. 29 Wei, J., 2000. Point-condensations generated by the Gierer-Meinhardt system: a brief survey in “New trends in nonlinear partial differential equations”. (Y Morita, H Ninomiya, E Yanagida and S Yotsunani editors), pp. 46–59. Wu, J., Berland, K. M., 2008. Propagators and time-dependent diffusion coefficients for anomalous diffusion. Biophys. J. 95(4), 2049–2052. Xiao, F., Hrabe, J., Hrabetova, S., 2015. Anomalous extracellular diffusion in rat cerebellum. Biophys. J. 108, 2384–2395. Zhang, Y., Green, C. T., Tick, G. R., 2015. Peclet number as affected by molecular diffusion controls transient anomalous transport in alluvial aquifer-aquitard complexes. J. Contaminant Hydrology 177-178, 220–238. 30 D D1 D2 τfin 0 τ 1 spike pattern 2 spike pattern (0) H0 (0) (0) (0) (0) H0 (τ ) H0 (τ ) H0 (0) (0) H0 (2τ× ) (0) H0 (2τ× ) x0 (0) x0 (τ ) x0 (2τ× ) x0 (0) x0 (τ ) x0 (2τ× ) Figure 3: Piecewise constant approximation to a smooth transition. Upper panel: D(τ ) based on (16b), dashed curve, and (0) its step-discontinuous discretisation in ten equal steps (solid lines). Bottom panels: phase plane trajectory H0 (x0 ) for one (left) and two (right) spike constellations. Line styles conform to curves in upper panel. Chemical exponents used (p, q, m, s) = (2, 1, 2, 0), with qualitatively similar results for sets (2,1,3,0), (3,2,2,0) and (4,2,2,0). Spike width ǫ = 0.1. 31 (0) H0 eq,D1 (0) eq,D2 Hi (0) H0 x0 (τi ) x0eq x0 (τ× ) x0 D D1 D2 0 τ× τ (0) Figure 4: Three spike pattern transition for varying diffusivity. Upper panel: spike height H0 (τ ) (outer spike, starting at circle (0) marks H0 (τi ) Di (0) (0) , i={1,2}) and H1 (τ ) (central spike, starting at diamond marks H1 (τi ) Di , i={1,2}) versus locus position x0 (τ ) for inner layer (solid) and step discontinuous (dashed) functions (equations (16b) with ε = 2 and (16a) respectively). Dotted lines show drift towards equilibria (cross marks) for corresponding constant diffusivities D1 and D2 . Lower panel: time dependence D τ . Chemical exponents used (p, q, m, s) = (2, 1, 2, 0), with qualitatively similar results for sets (3,2,2,0) and (0) (0) (4,2,2,0), for the latter with a minor difference of H1 (τi ) D1 > H0 (τi ) D2 . Spike width ǫ = 0.1. (0) H0 eq,D1 (0) eq,D2 Hi (0) H0 x0 (τi ) x0eq x0 (τ× ) x0 D D1 D2 0 τ× τ (0) Figure 5: An example of non-monotonicity for a three spike pattern transition. Upper panel: spike height H0 (τ ) (outer spike, (0) starting at circle marks H0 (τi ) (0) Di (0) , i={1,2}) and H1 (τ ) (central spike, starting at diamond marks H1 (τi ) Di , i={1,2}) versus locus position x0 (τ ) for inner layer (solid) and step discontinuous (dashed) functions (equations (16b) with ε = 2 and (16a) respectively). Dotted lines show drift towards equilibria (cross marks) for corresponding constant diffusivities D1 and  D2 . Lower panel: time dependence D τ . Chemical exponents used (p, q, m, s) = (2, 1, 3, 0). Spike width ǫ = 0.1. 32 H10 1 H10 (τi ) H10 (τi ) D1 D2 3 x0 x1 max τ× x0 x1 max . D1 x0 x1 x0 x1 D2 Figure 6: Four spike pattern transition for varying diffusivity: inner to outer spike height ratio H10 (τ ) versus locus position  ratio x0 (τ ) x1 (τ ) for inner layer (solid) and step discontinuous (dashed) functions (equations (16b) with ε = 2 and (16a) respectively). Dotted lines show drift towards equilibrium (cross mark) for corresponding constant diffusivities D1 and D2 . Time dependence D τ same as in figure 4. Chemical exponents used (p, q, m, s) = (2, 1, 2, 0), with qualitatively similar results for sets (2,1,3,0), (3,2,2,0) and (4,2,2,0). Spike width ǫ = 0.1. (0) H0 (τ ) 0.25 0.19 −0.5 0 x0 (τ ) (0) Figure 7: Single spike transition for varying diffusivity with sub-diffusion: spike height H0 (τ ) versus locus position x0 (τ ) for inner layer (solid) and step discontinuous (dashed) functions (equations (16b) with ε = 2 and (16a) respectively). Anomaly exponents: γ = 1 (normal, thick black lines), γ = 0.25 (thin black lines, slowest transition), γ = 0.15 (thick grey lines) and γ = 0.1 (thin grey lines). Chemical exponents used (p, q, m, s) = (2, 1, 2, 0), with qualitatively similar results for sets (2,1,3,0), (3,2,2,0) and (4,2,2,0). Spike width ǫ = 0.1. Time dependence D τ same as in figure 1. 33 (0) H0 (τ ) 0.25 0.19 −0.5 0 x0 (τ ) (0) H0 (τ ) 0.25 0.19 −0.5 0 x0 (τ ) Figure 8: Comparison of full numerical (solid) and asymptotic (dotted) solutions for a single spike transition for varying diffu(0) sivity with sub-diffusion: spike height H0 (τ ) versus locus position x0 (τ ) for inner layer (upper panel) and step discontinuous (lower panel) functions (equations (16b) with ε = 2 and (16a) respectively). Anomaly exponents: γ = 1 (normal, thick black lines), γ = 0.25 (thin black lines, maximal error during transition), γ = 0.15 (thick grey lines) and γ = 0.1 (thin grey lines). Chemical exponents used (p, q, m, s) = (2, 1,  2, 0), with qualitatively similar results for sets (2,1,3,0), (3,2,2,0) and (4,2,2,0). Spike width ǫ = 0.1. Time dependence D τ same as in figure 7. (0) H0 (τ ) 0.153 0.136 −0.75 −0.5 x0 (τ ) (0) Figure 9: Two spike pattern transition for varying diffusivity with sub-diffusion: spike height H0 (τ ) versus locus position x0 (τ ) for inner layer (solid) and step discontinuous (dashed) functions (equations (16b) with ε = 2 and (16a) respectively). Anomaly exponents: γ = 1 (normal, thick black lines), γ = 0.25 (thin black lines), γ = 0.15 (thick grey lines) and γ = 0.1 (thin grey lines). Chemical exponents used (p, q, m, s) = (2, 1, 2, 0), with qualitatively similar results for sets (2,1,3,0), (3,2,2,0) and (4,2,2,0). Spike width ǫ = 0.1. Time dependence D τ same as in figure 2. 34 H10 (τ ) 1 0.88 -0.8 x0 (τ ) − 32 Figure 10: Three spike pattern transition for varying diffusivity with sub-diffusion: inner to outer spike height ratio H10 (τ ) versus locus position x0 (τ ) for inner layer (solid) and step discontinuous (dashed) functions (equations (16b) with ε = 2 and (16a) respectively). Anomaly exponents: γ = 1 (normal, thick black lines), γ = 0.25 (accelerated transition, thin black lines, leftmost pair of curves), γ = 0.15 (delayed transition, thick grey lines, dashed line indistinguishable from γ = 1 case) and γ = 0.1 (extremely delayed transition, thin grey lines, rightmost pair of curves). Chemical exponents used (p, q, m, s) = (2, 1, 2, 0), with qualitatively similar results for sets (2,1,3,0) and (4,2,2,0). Spike width ǫ = 0.1. Time dependence D τ same as in figure 9. Set (3,2,2,0) gave results similar to the two spike case, not shown. . x0 x1 max xx10 max x0eq + δ0 x1eq + δ1 3 0 1.5 τ Figure 11: Four spike pattern transition  for varying diffusivity with γ = 1: comparison of full system (63) (solid) and linear system (64) (dashed) drift solution x0 x1 versus τ for inner layer (black) and step discontinuous (grey) functions (equations (16b) with ε = 2 and (16a) respectively). Dotted curves show drift with no transition, corresponding to D1 and D2 . Chemical exponents used (p, q, m, s) = (2, 1, 2, 0), with qualitatively similar results for sets (2,1,3,0), (3,2,2,0) and (4,2,2,0). Spike width ǫ = 0.1. 35 . x0 x1 3 0 τ× τ Figure 12: Four spike pattern transition  for varying diffusivity: comparison of full system (63) (solid) and partly linearised system (64) (dashed) drift solution x0 x1 versus τ for inner layer (black) and step discontinuous (grey) functions (equations (16b) with ε = 2 and (16a) respectively) and two values of anomaly exponent: γ = 0.95 (thick curves) and γ = 0.9 (thin curves). Chemical exponents used (p, q, m, s) = (2, 1, 2, 0), with qualitatively similar results for sets (2,1,3,0), (3,2,2,0) and (4,2,2,0). Spike width ǫ = 0.1. H10 1 3  x0 x1 Figure 13: Four spike pattern transition for  varying diffusivity and anomaly exponent γ < 1 : inner to outer spike height ratio H10 (τ ) versus locus position ratio x0 (τ ) x1 (τ ) for inner layer (solid) and step discontinuous (dashed) functions (equations (16b) with ε = 2 and (16a) respectively): γ = 0.95 (thick black curves), γ = 0.85 (thin black curves), γ = 0.75 (thick grey curves), γ = 0.65 (thin grey curves). Chemical exponents used (p, q, m, s) = (2, 1, 2, 0), with qualitatively similar results for sets (2,1,3,0), (3,2,2,0) and (4,2,2,0). Spike width ǫ = 0.1. 36