Math. Model. Nat. Phenom. Vol. #, No. #, 2017, pp. 0–29 DOI: 10.1051/mmnp/2017*** Transitions in asymmetric spike patterns with anomalous diffusion Y. Nec1 ∗ 1 Thompson Rivers University, British Columbia, Canada Abstract. Complex natural systems at times manifest transitions between disparate diffusive regimes. Efforts to devise measurement techniques capable of identifying the cross-over moments have recently borne fruit, however interpretation of findings remains contentious when the bigger picture is considered. This study uses a phenomenological reaction – diffusion model to obtain a sufficiently general, yet tractable description for the dynamics of a pattern qualitatively redolent of molecular clusters subject to transient anomalous diffusion mechanisms. The resulting system of equations is investigated for substantiation of the reasons difficulties are encountered when attempting to distinguish between various diffusive regimes in experimental settings: a non-monotonic dependence of the pattern’s evolution on parameters defining the diffusion mechanism is a common occurrence, as is a non-injective mapping between a given sequence of diffusion regimes and ensuing drift behaviour. Albeit formally the developed description is limited to piecewise constant anomaly functions, it is useful in computation of approximations to smooth realistic transitions. Keywords and phrases: fractional differential equation, variable order differential equation, matched asymptotic expansions, non-linear pattern formation, reaction – diffusion, Gierer-Meinhardt model, spike solution, transient diffusion regimes, anomalous diffusion Mathematics Subject Classification: 34E10, 34E20, 35K57, 35K67, 35K92, 35Q92, 35R11, 37M05 1. Introduction Anomalous diffusion is a collective name for dispersion processes wherein no linear growth of the mean square displacement is observed, i.e. ⟨r2 ⟩ ≁ t [Klages et al., 2008]. One well explored type of anomaly manifests ⟨r2 ⟩ ∼ tγ with the index γ ranging 0 < γ < 1, the sub-linear growth having granted it the name of sub-diffusion. Theoretically it can be obtained by normal diffusion over a fractal domain [ben Avraham and Havlin, 2000] or over a domain of integer dimensionality via a memory operator such as a fractional time derivative [Metzler and Klafter, 2000]. In experiments this dependence with a constant index was often observed over limited periods of time, thereafter the index indicating transition to a distinct anomalous regime or reverting to normal diffusion, γ = 1. For instance, Saxton [1994] discusses anomalous to normal cross-over due to the crowdedness of such environments as living cell cytoplasm and organelles. A later work by the same author further opines that in the absence of thermal equilibration the ensuing regime will remain anomalous with a different index [Saxton, 2007]. A continuous variation of the anomaly index with crowdedness is shown in Banks and Fradin [2005]. Specific values to be expected in such systems were measured in Sanabria et al. [2007]. Dentz et al. [2004] and Bijeljic et al. [2011] explored flow through a heterogeneous porous medium, another system, where such transitions might be expected, theoretically corresponding to a change in the dimension of the underlying fractal domain. One more interesting example is that of a single file diffusion, where the anomaly index can be obtained analytically as γ = 12 [Hahn and Kärger, 1995]. As remarked therein, if the domain is such that the one dimensional constraint holds over only some portion thereof or is altogether removed at a given point in time, regular diffusion is obtained, begetting a cross-over. Lateral diffusion within a membrane or generally an envelope, preceded or followed by motion in the cytoplasm or other less confined compartment, conceivably conforms to such a scenario and was investigated in Cabal et al. [2006] as well as Bakalis et al. [2015]. Inspired by the foregoing set of works, the current contribution aims to investigate transitions between diffusion regimes in the Gierer-Meinhardt model, a phenomenological reaction – diffusion system known to admit solutions in the form of patterns of spikes both with regular diffusion and sub-diffusion. In essence a spike pattern is a solution comprising a set of high concentration pulses, stationary or drifting. A conceptual analogy can be drawn between these localised high concentrations of reactants and clusters of ligands moving and reacting within a cell, thereby permitting to explore the possible transient dynamics between disparate diffusion regimes without the onus of modelling specific intricate molecular processes. With regular diffusion variegated spike constellations were extensively studied, cf. Iron et al. [2001], Iron and ∗ Corresponding author. E-mail: cranberryana@gmail.com © EDP Sciences, 2017 Y. Nec Transitions in anomalous spike patterns Ward [2002] and references therein. With sub-diffusion the existence of spike solutions was established in Nec and Ward [2012]. More recently transient behaviour was studied in sub-diffusive patterns with a fixed diffusion coefficient [Nec, 2018], as well as a fixed anomaly exponent and reflective symmetry only [Nec, 2017]. These works established the feasibility of relaxing the constraint of a constant temporal derivative order or constant diffusion coefficient respectively, but not both. The existence of a pattern of spikes was explored and some limitations were discovered. The main focus was on attainment of closed form solutions for simplified constellations. The resulting systems were used for basic observations of behaviour within a limited parameter space, however no attempt was made to draw deeper physical inferences. Here combined transitions are explored for fully asymmetric patterns, elucidating the similarities and differences between the system responses to changes in diffusion coefficient and anomaly index. Amongst other comprehensions, these findings shed light on the general difficulty of classifying experimental observations as anomalous sub-diffusion versus slowed down normal diffusion in short term realistic processes [Soula et al., 2013]. In §2 the Gierer-Meinhardt model is generalised to allow transitions in the anomaly index as well as diffusion coefficient, a reduced algebraic-differential system is derived for a constellation of spikes at quasi-equilibrium, formally valid for slowly varying diffusivity function of unity order of magnitude and piecewise constant anomaly index. Transition scenarios are discussed in §3. Numerical solutions for instantaneous cross-over events are presented in §4, followed by construction of approximations for realistic smooth transitions. §5 concludes with a summary of the findings. Spike solutions are unstable under certain parameter regimes. The predominant parameters that affect the onset of instability are diffusivity D (amplitude stability, consult Iron et al. [2001], Iron and Ward [2002], and Ward and Wei [2002] for detail) and reaction time constant τo (oscillatory, competition as well as dynamically triggered instability [Ward and Wei, 2005],[Ward and Wei, 2003],[Tzou et al., 2015]), but also kinetic exponents and the number of spikes in a constellation. The choice of parameters for all numerical solutions presented hereinafter is well within the boundaries of known stability domains. Comparison of numerical solutions to the fractional partial differential system and the reduced algebraic-differential system remains an open problem. So are stability domains within the parameter space and modes of instability onset subject to variation in D and γ. Of particular interest is the possibility of instability triggered dynamically in the course of a transition. 2. Generalised Gierer-Meinhardt model The generalised Gierer-Meinhardt model with Neumann boundary conditions to be studied reads ∂tγ a = ǫ2 axx − a + ap hq − 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)hxx − h + ǫ−1 ax (±1, t) = hx (±1, t) = 0 , (2.1a) (2.1b) (2.1c) wherein a(x, t) and h(x, t) denote activator and inhibitor concentrations respectively. The inhibitor diffusivity D varies slowly in time, whilst always remaining of unity order of magnitude, i.e. D ∼ O(1) ∀ t. If D ∈ C 1 , the slow time d variation implies D ∼ O(ǫα ). The scale α > 0 is related to γ and discussed hereunder. The activator diffusivity involves dt an asymptotically small parameter 0 < ǫ ≪ 1. The reaction time constant τo > 0, the quadruple of kinetic exponents (p, q, m, s) satisfies p−1 m p > 1, q > 0, m > 0, s ⩾ 0, < , (2.1d) q s+1 and the anomaly index γ ranges 0 < γ ⩽ 1. When γ = 1, system (2.1) comprises two conventional partial differential equations. When 0 < γ < 1, the derivative of order γ is defined (for a univariate function) as 1 dγ f (t) = − dtγ Γ (−γ) ∫0 t f (t) − f (t − ζ) dζ, ζ γ+1 0 < γ < 1, (2.1e) where Γ denotes the Gamma function. The integration being with respect to ζ, the order γ is easily generalised to be a time dependent function γ(t), whereby (2.1e) becomes a variable order memory operator. Multifarious operators of that ilk exist, at times incorporating differentiation with respect to t or imbedding the delay through γ(t − ζ). In those cases relaxation of the constant order constraint requires more caution, cf. Ramirez and Coimbra [2009]. Infra the solution is constructed for a piecewise constant function γ(t) = ∑ γj {θ(t − t× j−1 ) − θ(t − t× j )}, j j = 1, 2, . . . , (2.2) with 0 < γj ⩽ 1 constant anomaly indices, θ(t) the Heaviside step function, t×j corresponding cross-over times satisfying t×j > t×j−1 and t×0 = 0 without loss of generality. If γ(t) ≡ 1, the diffusion mechanism is the regular Fickian diffusion with linear dispersion ⟨r2 ⟩ ∼ t, recovering the classical model [Gierer and Meinhardt, 1972]. With 0 < γ(t) < 1 the dispersion is sub-diffusive, i.e. ⟨r2 ⟩ ∼ tγ(t) . An arbitrary number of transitions is possible, and it is assumed that t×j − t×j−1 ∼ O(ǫ−α ), to wit γ varies slowly on the same time scale as D. A plain process is obtained with all γj equal. The choice γj = 1, 1 Y. Nec Transitions in anomalous spike patterns 0 < γj+1 < 1 gives a regular to anomalous cross-over and vice versa, whilst 0 < γj < 1, γj ≠ γj+1 yields transitions between distinct sub-diffusive regimes. Generally the cross-over points of γ(ǫα t) are not correlated with the functional form of D(ǫα t). Quondam studies used an equivalent form of system (2.1) with a different small parameter ǫ, whose scale involved the anomaly index γ [Nec and Ward, 2012], [Nec, 2018]. Identically to any asymptotic formalism there is a degree of freedom in the choice of power for ǫ. Due to the memory operator (2.1e) at least one out of the spatial and temporal asymptotic scales must be a function of γ. The current scalings are chosen so that only the temporal scale depends thereon. 2.1. Spike solutions A spike pattern comprises a constellation of n pulses centred at arbitrary loci xi , {i = 0, . . . , n − 1} within the domain. The outer solution for the activator away from xi is the background state whereas for the inhibitor expand a(x, t) ≡ 0, (2.3a) h ∼ h(0) (x, ǫα t) + o(1). (2.3b) Near the spike locus the inner spatial variable must be taken as yi = (x − xi )/ǫ to make the diffusion term in (2.1a) of order O(1). The corresponding slow time variable is τ = ǫα(t) t or t = τ ǫ−α(τ ) , the inversion in this simple form possible due to γ(t) being piecewise constant (equation (2.2)). The slow time cross-over moments are τ×j = ǫα t×j , j = 1, 2, . . . Therefore henceforth γ and α are regarded as functions of t or τ in accord with context. The inner asymptotic solutions thus are def Ai (yi , τ ) = a(xi + ǫyi , τ ǫ−α(τ ) ) ∼ Ai (yi , τ ) + ǫAi (yi , τ ) + ⋯ (2.4a) Hi (yi , τ ) = h(xi + ǫyi , τ ǫ−α(τ ) ) ∼ Hi (yi , τ ) + ǫHi (yi , τ ) + ⋯ . (2.4b) (1) (0) (0) (1) Expansions (2.4) are somewhat unorthodox due to the time dependence of α, nonetheless valid for any step-like function as in (2.2): at any cross-over moment t× j the scales of the entire asymptotic series are instantaneously adjusted to match the order of the memory operator. Hereafter the course of solution follows closely the construction of a spike pattern formerly devised for constant diffusivity D [Iron and Ward, 2002] and anomaly index γ [Nec and Ward, 2012], with the purport to unify the method for a piecewise constant γ(ǫα t) and generic D(ǫα t) ∼ O(1). Substituting (2.4) into (2.1) and collecting O(1) terms, (0) p Ai (0) (0) 2 ∂yi Ai − Ai + (0) q = 0, Hi D(τ ) ∂y2i Hi (0) = 0, = 0, (2.5a) = lim ± h(0) . (2.5b) (0) lim Ai ∣yi ∣Ð→∞ (0) lim Hi yi Ð→±∞ xÐ→xi Note that from D(ǫα t) ∼ O(1) it follows that D(τ ) ∼ O(1) ∀ τ , whence no inconsistence ensues when forbearing to expand D(τ ) in apposite powers of ǫ. Hi must be bounded as ∣yi ∣ Ð→ ∞, thus equation (2.5b) implies that Hi a function of τ . Therefore (2.5a) is solved as (0) (0) Ai = Hi u(y) = ⎛p + 1⎞ ⎝ 2 ⎠ (0) wherein the function u is the homoclinic solution to (0) β (τ )u(yi ), β= 1/(p−1) u′′ − u + up = 0, cosh−2/(p−1) must be solely q , p−1 (2.6a) ⎛p − 1 ⎞ y ⎠ ⎝ 2 (2.6b) lim u(y) = 0, ∣y∣Ð→∞ (2.6c) decaying exponentially and thus matching correctly the outer solution (2.3a). Under the premise of quasi-equilibrium the spike drifts slowly, to wit the locus xi depends on τ yi (t) ≡ (x − xi (τ ))/ǫ. (2.7) Thus at the next order the temporal derivative on the left hand side of (2.1a) contributes a term. The scale α must be chosen so that said term balances the correction terms on the right hand side of (2.1a). Since γ is piecewise constant, away from the cross-over moments τ×j the dependence γ(τ ) does not affect the balance, whilst at τ×j the scales are instantaneously (0) adjusted. With γ = 1 (regular diffusion) differentiation of Ai (yi (τ ), τ ) with respect to the first argument and invocation dxi (0) . The derivative with respect to the second argument of Ai (yi (τ ), τ ) is of a of the chain rule yield the drift velocity dτ 2 Y. Nec Transitions in anomalous spike patterns smaller order of magnitude. For γ < 1 (sub-diffusion) a similar separation of scales holds, only no elegant generalisation of the chain rule exists for operator (2.1e), therefore a specific asymptotic result was derived (consult Appendix A): dxi dxi γ−1 γ ∣ ∣ Dyi Ai (yi ) , dτ dτ ∂τγ Ai (yi (τ )) ∼ −ǫ−γ (2.8a) ⎧ ⎫ ⎪ ⎛ dxi ⎞⎪ ⎪ ⎪ dξ ⎨Ai (yi ) − Ai yi + ξ sgn ⎬ γ+1 . (2.8b) ⎪ dτ ⎠⎪ ξ ⎝ ⎪ ⎪ 0 ⎩ ⎭ This result was formerly obtained for a constant anomaly index γ and different, but equivalent scalings [Nec and Ward, 2012]. Here it is confirmed to hold for the more general case of a piecewise constant γ. Comparison of the order of magnitude of correction terms on the right hand side of (2.1a), O(ǫ), with αγ − γ on the left, as ensues from combining 1 dxi Dγyi Ai (yi ) ≡ sgn Γ (−γ) dτ ∞ ∫ ∂tγ f (t) = ǫαγ ∂τγ F (τ ) ∀ γ(t), F (τ ) ≡ f (ǫ−α τ ) (2.9) with (2.8a), yields the balance αγ − γ = 1, or α = 1 + 1/γ. With regular diffusion (γ = 1) the diffusivity ratio ǫ2 and slow time scale ǫ1+1/γ coincidentally equal. Substituting (2.4) into (2.1), using (2.8) and collecting O(ǫ) terms for the activator equation (2.1a) and O(ǫ−1 ) for the inhibitor equation (2.1b), one obtains q dxi dxi γ−1 γ (1) p−1 ⎞ (1) p ∣ ∣ Dy i u , − 1 + pu A = u H − i (0) dτ dτ ⎠ i ⎝ ∂yi2 Hi (2.10a) ∂ 2 (1) (0) βm−s m u . Hi = −Hi 2 ∂yi (2.10b) (0) −β ⎛ ∂ Hi 2 D(τ ) Since in this setting 0 < γ ⩽ 1, for 0 < γ < 1 the operator Dγyi u is defined by the fractional integral in (2.8b), whereas for ≡ u′ (yi ). The term τo ∂tγ Hi γ = 1 it is taken as Dγyi u∣ (0) (0) γ=1 in (2.1b) is negligible at this order since τo ∂tγ Hi ∼ O(ǫ1+γ ) if τo ∼ O(1) by (2.9) with α = 1 + 1/γ, or τo ∂tγ Hi ∼ o(ǫ−1 ), to hold whilst τo ∼ o(ǫ−(2+γ) ). By the Fredholm alternative, (2.10a) is soluble if ⎫ ⎧ ⎪ q du ⎪ dxi dxi γ−1 γ ⎪ ⎪ (1) (2.11a) ⎨ (0) up Hi − ∣ ∣ Dyi u⎬ dyi = 0. ⎪ ⎪ dyi ⎪ H dτ dτ ⎪ R ⎭ ⎩ i The first term in (2.11a) is integrated by parts twice and simplified with the aid of (2.6c) along with the facts that u and (1) ∂y2i Hi are even functions in yi by (2.6b) and (2.10b) respectively. Thence (2.11a) becomes (0) ∫ q (0) 2(p + 1)Hi ∫ R ∫ dy D u dy . ∂Hi (1) dxi dxi γ−1 ∂Hi (1) + lim )= − ∣ ∣ yi Ð→−∞ ∂yi yi Ð→∞ ∂yi dτ dτ up+1 dyi ( lim du R i γ yi i (2.11b) To obtain the leading order equation for the outer solution substitute (2.3b) into (2.1b) and collect terms of O(1): n−1 D(τ ) ∂x2 h(0) − h(0) = −bm ∑ Hi (0) βm−s i=0 δ(x − xi ), bm = ∫ um dy, R (2.12a) wherein the weighted sum of Dirac delta functions δ(x − xi ) on the right hand side appears due to the localised nature of the non-linear reaction term: 1 am n−1 (2.12b) ∼ ∑ bi δ(x − xi ) ǫ hs i=0 with 1 bi = ǫ ∫ x+ i am dx = hs − xi ∫ ∞ −∞ am (0) βm−s dyi ∼ Hi s h ∫ u dy . m (2.12c) R Equation (2.12a) is solved as h(0) (x, τ ) = ∑ ci (τ )G(x; xi ; τ ), n−1 i=0 ci = −bm Hi (0) βm−s (2.13a) with G(x; xi ; τ ) being the Green’s function (prime denoting differentiation with respect to the first argument) D(τ ) G′′ − G = δ(x − xi ), G′ (±1; xi ; τ ) = 0, (2.13b) somewhat unconventional due to the diffusivity D(τ ) depending on the slow time scale τ . The computation of the left hand side in (2.11b) is completed with (1) lim yi Ð→−∞ ∂Hi ∂yi + lim yi Ð→∞ (1) ∂Hi ∂yi = lim − ∂x h(0) + lim + ∂x h(0) = xÐ→xi 3 xÐ→xi (2.14) Y. Nec Transitions in anomalous spike patterns ⎫ ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎞ 1 (0) βm−s ⎛ ′ − ⎪ ⎪n−1 (0) βm−s ′ ′ + G (xi ; xi ; τ ) + G (xi ; xi ; τ ) ⎬ . G (xi ; xj ; τ ) + Hi −2bm ⎨ ∑ Hj ⎪ 2 ⎝ ⎠⎪ ⎪ ⎪ j=0 ⎪ ⎪ ⎪ ⎪ ⎭ ⎩ j≠i Further matching of leading order inner and outer inhibitor solutions h(0) (xi , τ ) = Hi (τ ) (0) (2.15) in (2.13a) yields the algebraic-differential system for the spikes’ heights and drift of loci: ⎧ ⎫ ⎪ ⎪ ⎪ n−1 ⎪ ⎪ βm−s dxi dxi γ−1 qbm f (p; γ) ⎪ ⎪1 (0) βm−s ′ − ⎪ (0) ′ ′ + ⎨ G (x ; x ; τ ) ⎬, H (G (x ; x ; τ ) + G (x ; x ; τ ))+ ∣ ∣ = H ∑ i j i i i i j i (0) ⎪2 ⎪ dτ dτ ⎪ (p + 1)Hi ⎪ j=0 ⎪ ⎪ ⎪ ⎪ j≠i ⎩ ⎭ n−1 Hi (τ ) = −bm ∑ Hj (0) f (p; γ) ≡ ∫ u p+1 G(xi ; xj ; τ ) , (2.16b) u′ (y)Dγy u dy (2.16c) (0) βm−s j=0 dy R /∫ R (2.16a) with f (p; γ(τ )) varying in accord with γ. When γ = 1, f (p; 1) = 2(p + 1)/(p − 1). The functional dependence of f (p; γ) is given in figure 1. During a transition in γ an instantaneous adjustment f (p; γj ) ↝ f (p; γj+1 ) will occur. f (p; γ) 9 p=2 6 p=3 4 p=5 p=4 3 0 γ 1 Figure 1. Factor f (p; γ) in algebraic-differential system (2.16). Cross marks show the normal values f (p; 1) = 2(p+1)/(p−1). f (p; γ) Ð→ ∞ as γ Ð→ 0+ . Reproduced from figure 3 of Nec [2018] for convenience. It remains to solve (2.13b) as √ ⎧ ⎪ ⎪ 1 cosh ((xi − 1)/ D ) x+1 ⎪ ⎪ ⎪ −√ cosh √ ⎪ √ ⎪ ⎪ D D ⎪ sinh (2/ D ) ⎪ ⎪ ⎪ ⎪ G(x; xi ; τ ) = ⎨ √ ⎪ ⎪ ⎪ ⎪ cosh ((xi + 1)/ D ) ⎪ 1 x−1 ⎪ ⎪ ⎪ −√ cosh √ ⎪ √ ⎪ ⎪ D D ⎪ sinh (2/ D ) ⎪ ⎩ x < xi (2.16d) x > xi with D varying in τ if not taken constant. In summary, system (2.16) is the reduction of the fractional partial differential system (2.1) for the drift of a stable constellation of n spikes. The reduction is valid for a piecewise constant function 0 < γ(t) ⩽ 1 with any number of 4 Y. Nec Transitions in anomalous spike patterns transitions occurring at O(ǫ−(1+1/γ) ) intervals of the time t and D(τ ) ∼ O(1). At each discontinuity point of γ the drift time scale is adjusted. System (2.16) is unsuitable for tracking evolution near and past the onset of instability: D(τ ) must remain well within the bounds of the stability domain at all times. Albeit the foregoing formalism is valid only for piecewise constant γ(τ ), it is shown hereinafter that (2.16) is viable in approximation of the solution for any function γ(τ ) amenable to a piecewise constant discretisation. D(τ ) may have jump discontinuities. The pattern state immediately prior to the discontinuity is to be regarded as the initial condition thereupon. If D is smooth, infinitesimal relaxation times in response to the respective increments in D eventuate. As an example the drift of a single spike was simulated via system (2.1) with D(τ ) given by (3.3a) for γ ≡ 1, and compared to the solution as given by the reduced system (2.16). Figure 2(a) depicts the ensuing evolution. Some explanations are in order. The parameters were chosen so as to illustrate the source and magnitude of the error involved. The simulation begins at t = 0 with the spike height corresponding to the initial position as predicted by (2.16b). It is imperative that the initial state be the correct quasi-equilibrium solution, as otherwise the relaxation would involve extraneous effects. Provision of a correct quasi-equilibrium height is possible for one spike only due to the non-linearity of (2.16b). At each moment in time the spike location is extracted by arg max a(x, t). Equation (2.6a) is then used to find x (0) H0 : since the homoclinic (2.6b) is neither time nor locus dependent, max a(x, t) ∼ H0 (0) β x 1/(p−1) (0) β ⎛ p + 1 ⎞ max u(y) = H0 y ⎝ 2 ⎠ (2.17a) is expected to hold throughout the simulation, yielding −1/(p−1) 1/β ⎛ ⎞ ⎛p + 1⎞ (0) ⎟ . H0 (t) ∼ ⎜max a(x, t) x ⎝ 2 ⎠ ⎝ ⎠ (2.17b) The correspondence between the solid and dashed curves is obvious, to wit both systems faithfully capture the transition. Moreover, the computation for (2.16) with the curve (3.3a) discretised into 200 equal steps yielded a virtually indistinguishable curve. The visible deviation stems from the finiteness of ǫ: the initial state corresponds to the limit ǫ Ð→ 0, therefore the first phenomenon observed is relaxation to the height conforming to the finite value of ǫ used. This slight difference in height is thence maintained throughout the evolution. This error diminishes with ǫ, as is depicted by figure 2(b). Note that the error bound for the height is O(ǫ) by (2.4b). In fact the reduced system gives a better accuracy than 1 expected, since the error never exceeds 10 ǫ. Mere shifting of the full solution curve to compensate for the initial height error and relaxation time results in a phase path nearly overlapping the prediction of (2.16). In light of the above (2.16) reconstructs the qualitative behaviour of the pattern with a fair accuracy, better than is expected from the formal error bound of the asymptotic solution (albeit this might deteriorate with γ < 1). It stands to reason there would be aspects of the dynamics of (2.1) that the algebraic-differential system would be unable to capture: there will always be an O(1) time interval after a jump discontinuity in D(t) over which systems (2.1) and (2.16) do not agree. In particular consecutive discontinuities must be separated by a relaxation period of O(ǫ−(1+1/γ) ) in t. System (2.16) might be able to anticipate response delay upon instantaneous transitions for constellations with a sufficient number of spikes, cf. figure 6 of Nec [2017]. Not so for a single spike or equivalent (e.g. the dimensionality of a symmetric pair is equal to a single spike). Figure 3 sets forth a classic example thereof. Unsurprisingly the full solution manifests discernible relaxation time, whilst the reduced system predicts an instant cross-over, the leading order approximation entailing an adjustment to the required quasi-equilibrium state. However, bearing in mind the foregoing remarks on the prevised error bound, the difference is minor. Inspection of (2.16) permits several qualitative inferences concerning the response to variation in D and γ. Changes in D affect the spike height and velocity by (2.16b) and (2.16a) respectively, whereas a shift in γ would have an immediate impact on the locus velocity alone. However, the height will respond by the coupling to (2.16b). The adjustment of velocity and hence height can be swift and easily observable if the jump in γ is sufficient – as is seen from figure 1, the factor f (p; γ) is a rapidly descending function. Therefore given the eventuating drift pattern only, it might be virtually impossible to predict whether the response is due to a transition in D or γ. 3. Transition scenarios System (2.16) is to be integrated subject to a piecewise constant anomaly index γ (consult equation (2.2)) in combination with an apposite variation of the diffusion coefficient D. The particular arrangement of the cross-over moments and relative changes of the two parameters is referred to as a drift scenario, and constitutes the salient focus of this contribution. Past experience indicates that abrupt, yet smooth front-like transitions in D have only a quantitative impact in the form of a slight delay near the cross-over moment, whereas step discontinuous D is the simplest form to integrate whilst successfully reconstructing the dynamics [Nec, 2017]. For this reason piecewise constant D(τ ) functions are implemented here. In a general physical system a shift in the diffusion coefficient D might stem from a molecular change that does not affect the species identity, i.e. within the model framework the reactant functions as the same species, e.g. as a result of interaction with a catalyst or alteration of its spatial configuration. By contrast an abrupt transition in the 5 Y. Nec Transitions in anomalous spike patterns anomaly index γ might be effectuated by an essential change in the medium, for instance moving into or leaving a more crowded compartment, i.e. one of a disparate effective fractional dimension. Whatsoever the underpinning physics, these are independent events, whence variegated sequences of transitions are of interest. Two scenarios were found to be the building blocks of any commonly observed dynamics: a single cross-over and a double cross-over such that the initial value of the parameter at hand is recovered. Obviously, the latter is a special case of a pair of single events, and any complex transition sequence can be composed of consecutive blocks of these types. Simultaneous transitions of both parameters contribute no qualitative novelty, as one of them is the predominant as (0) regards the system’s response. By (2.16) a discontinuity in spike heights Hi will ensue with a change in D, whereas a shift in γ will entail a discontinuity in the drift velocity. Upon a concurrent transition both will be instantaneously adjusted, with only one dominating the response, as depends on the state of the pattern at the cross-over moment. That much said, it is possible to tune the parameters to yield an equal parts or annihilation response to such a synchronous change. Since the purport of this study is to analyse generic phenomena, no search after specific initial conditions or transition sequences with rare outcomes was made. Each transition is always followed by sufficient evolution time, so as not to create a control problem, where it is no longer possible to match outcome to input. 3.1. Instantaneous events: permanent shift The first type of scenario treated is the simpler one, where γ and D shift abruptly once at distinct cross-over moments. The four scenarios tested are constructed so that a change in the first parameter is either augmented (shift in the same direction) or contrasted (shift in the opposite direction) by the second transition. Whilst the sequence content is of import, swapping the order made only for quantitative differences in the final state of the pattern. The following set of transitions was found to produce well quantifiable effects (see also figure 4). For reference hereunder each sequence is encoded with upward and downward arrows based on the direction of change for each parameter. t×0 = 0, γ = 1, γ↡ D↡ γ↡ D↟ γ↟ D↡ γ↟ D↟ t×1 = 15, 1 γ = , 4 ⎧ γ ⎪ ⎪  γ=⎨ ⎪ ⎪ ⎩ γ ⎧ γ ⎪ ⎪  γ=⎨ ⎪ ⎪ ⎩ γ ⎧ γ ⎪ ⎪  γ=⎨ ⎪ ⎪ ⎩ γ ⎧ γ ⎪ ⎪  γ=⎨ ⎪ ⎪ ⎩ γ t×2 = 30, D = 2, 0 < t < t×1 t > t×1 0 < t < t× 1 t > t×1 0 < t < t× 1 t > t×1 0 < t < t×1 t > t× 1 (3.1a) D = 1, ⎧ D ⎪ ⎪  D=⎨ ⎪ ⎪ ⎩ D ⎧ D ⎪ ⎪  D=⎨ ⎪ ⎪ ⎩ D ⎧ D ⎪ ⎪  D=⎨ ⎪ ⎪ ⎩ D ⎧ D ⎪ ⎪  D=⎨ ⎪ ⎪ ⎩ D 0 < t < t×2 t > t×2 0 < t < t× 2 t > t×2 0 < t < t× 2 t > t×2 0 < t < t×2 t > t× 2 (3.1b) (3.1c) (3.1d) (3.1e) 3.2. Instantaneous events: recovery sequences In this type of transition out of γ and D the parameter undergoing a step change at t× 1 recovers its initial value at t× 2 (denoted by upward and downward harpoons), subsequently the other parameter is permitted to change at t×3 (denoted by upward or downward arrow). The following set of transitions was found to produce well quantifiable effects (see also figure 5). (3.2a) t× 0 = 0, t× 1 = 15, t× 2 = 30, t× 3 = 45, γ = 1, 1 γ = , 4 D = 2, D ⥯ γ↟ ⎧ γ ⎪ ⎪  γ=⎨ ⎪ ⎪ ⎩ γ 0 < t < t× 3 D ⥯ γ↡ ⎧ γ ⎪ ⎪  γ=⎨ ⎪ ⎪ ⎩ γ 0 < t < t×3 t > t×3 t > t×3 6 D = 1, ⎧ D ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ D = ⎨ D ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ D ⎧ D ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ D = ⎨ D ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ D 0 < t < t×1 t× 1 < t < t× 2 (3.2b) t > t×2 0 < t < t×1 t× 1 < t < t× 2 t > t×2 (3.2c) Y. Nec Transitions in anomalous spike patterns γ ⥯ D↟ γ ⥯ D↡ ⎧ γ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ γ = ⎨ γ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ γ ⎧ γ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ γ = ⎨ γ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ γ 0 < t < t×1 t×1 < t < t×2 t > t× 2 0 < t < t×1 t× 1 < t < t× 2 t > t×2 ⎧ D ⎪ ⎪  D=⎨ ⎪ ⎪ ⎩ D 0 < t < t× 3 ⎧ D ⎪ ⎪  D=⎨ ⎪ ⎪ ⎩ D 0 < t < t×3 t > t×3 t > t×3 (3.2d) (3.2e) 3.3. Smooth transition An instantaneous transition should be regarded as a first approximation of a more realistic smooth cross-over. System (2.16) is formally valid for any function D(τ ), but only piecewise constant γ(τ ). Nonetheless the presented formalism is viable in solution of more realistic problems. With a sufficiently fine discretisation the aftermath is qualitatively as smooth as could be expected from continuous D(τ ) and γ(τ ). As an example the smooth transient functions below were considered: 1 1 τ − τ× D(τ ) = (D + D ) − (D − D ) tanh 1+1/γ , (3.3a) 2 2 ǫ εD 1 τ − τ× 1 γ(t) = (γ + γ ) − (γ − γ ) tanh 1+1/γ . 2 2 ǫ εγ (3.3b) The inner layer width parameters εD and εγ are independent of the spike width scale ǫ, and D⥯ , γ⥯ are constant as in (3.1a). The ranges D − D and γ − γ are divided into the desired number of steps that need not be equal, although that is the simplest variant. The specific cross-over forms (3.3) or discretisation do not encumber the genericity of the method. At each evolution stage created thereby system (2.16) is formally valid. Therefore the reduction of the variable order fractional differential system (2.1) to the algebraic-differential system (2.16) at the limit ǫ Ð→ 0 in fact provides a feasible approximation technique for any function γ(τ ) that is amenable to a representation by a sum of weighted step functions. The analytical justification is as follows. For initial input γ(τ×0 ), D(τ×0 ) and xi (τ×0 ), i = {0, . . . , n − 1}, the conforming spike heights are the solutions of (2.16b). The constellation is therefore in quasi-equilibrium, whence the system is allowed (0) to evolve over the time span of the first prescribed step-like stage (τ× 0 , τ× 1 ). The final state is {xi (τ× 1 ), Hi (τ× −1 )}. The time interval must satisfy τ× 1 − τ× 0 ∼ O(1). At τ× +1 the parameters γ and D attain new values. The height of all spikes and their velocity are adjusted by (2.16b) and (2.16a) respectively. The eventuating discontinuity depends on the increments D(τ× +1 ) − D(τ× −1 ) and γ(τ× +1 ) − γ(τ× −1 ), rendering the constellation in quasi-equilibrium again to evolve by (2.16). Figure 2 illustrated this process for γ ≡ 1 with (3.3a) and a finely discretised counterpart thereof. Coarsely discretised examples appear in figure 3 of Nec [2017]. Examples of more complicated transition scenarios are given in §4. 4. Numerical solution Equation (2.16a) was integrated numerically using the classical four step Runge-Kutta method, solving (2.16b) as required, with (2.16c) and (2.16d) adjusted in accord with the variation of γ and D respectively. Observe that (2.16a) is an ordinary dxi dxi γ−1 differential equation, admitting solution by regular numerical integration techniques: the term ∣ ∣ can be written dτ dτ γ dxi dxi as sgn ∣ ∣ , thence upon confirmation the right hand side bears the appropriate sign, a 1/γ root casts the equation dτ dτ dxi into the familiar form. The only technical difference is the necessity to keep track of sgn at all times. dτ The integration was performed for constellations of one to four spikes with parameters chosen so that stable motion ensued. Tracking evolution near onset of instability is beyond the scope of this study. To verify the code correctness quondam results for reflectively symmetric constellations reported in Nec [2017] were reproduced without dimensionality reduction due to symmetry. Results of Nec [2018] were reproduced with fixed D. Hereinafter transitions with asymmetric initial states are explored, being the generic choice of parameters responsible for outcomes unencumbered by artificial impositions. The dimensionality of the parameter space of system (2.16) is sufficiently high to make a complete presentation of results impracticable. For a fully asymmetric pattern of n spikes there are n loci and spike heights, whose evolution depends on chemical exponents (p, q, m, s), initial conditions and the chosen γ − D sequence, not the least part thereof being the cross-over moments. The solutions reported infra were found to be both of qualitative interest and feasible for visualisation. It is understood that some of the choices are arbitrary and other parts of the parameter space might be investigated in a dedicated numerical study. 7 Y. Nec Transitions in anomalous spike patterns 4.1. Instantaneous events: permanent shift System (2.16) was solved for one to four spike constellations, several common sets of chemical exponents, various asymmetric initial conditions and four transition sequences (figure 4). When the properties of the resulting drift dynamics were juxtaposed with those of the more complex recovery sequences (figure 5), it was observed that they were a subset thereof. For this reason only a basic behaviour of a single spike is presented in figure 6 by way of introduction of the trends to be expected. The following typical transient dynamics become evident. † An increase (diminution) of γ entails a tangible magnification (reduction) in the rate of change of the spike height, as is seen from panel 6(a) by comparison of the slope on two sides of t× 1 for sequences with γ ↡ (equations (3.1b) and (3.1c), thick curves) versus those with γ ↟ (equations (3.1d) and (3.1e), thin curves). The drift velocity responds respectively, observed from identical scrutiny of the curves’ slopes in panel 6(b). † An increase (diminution) of D brings about a jump discontinuity in height up (down), as is seen from panel 6(a) by comparison of the level on two sides of t×2 for sequences with D ↡ (equations (3.1b) and (3.1d), black curves) versus those with D ↟ (equations (3.1c) and (3.1e), grey curves). By contrast to the height evolution, the response of the drift velocity to changes in D is more complicated. Examination of panel 6(b) establishes that the velocity might follow any trend: respond conversely to the change in D, e.g. the spike speeds up due to a drop in D (scenario γ ↡ D ↡, equation (3.1b), thick black curve), or slows downs with a rise in D (scenario γ ↡ D ↟, equation (3.1c), thick grey curve); but also exhibit virtually no response, as happens for γ ↟ D ↡ (equation (3.1d), thin black curve) and γ ↟ D ↟ (equation (3.1e), thin grey curve). The responses due to variation in D might be regarded as normal in the sense that they could be predicted from (2.16) if solved with the germane constant values of D. The non-linearity of (2.16b), to wit the compound βm − s, in conjunction with the non-linear dependence of G(x; xi ; τ ) on D (equation (2.16d)) render the trends non-trivial to trace analytically. Thereon is superimposed the impact of sub-diffusion. Scrutiny of (2.16a) uncloaks two opposing influences. The factor f (p; γ) is monotonic in γ (bar a narrow region near γ = 1 for p small enough), increasing rapidly as γ diminishes and thus being instrumental in speeding the spike up. On the other hand, the magnitude of the right hand side of (2.16a) remains dxi dxi γ−1 dxi dxi γ below unity at all times, whence upon re-writing ∣ ∣ = sgn ∣ ∣ the 1/γ > 1 root curbs the foregoing effect. dτ dτ dτ dτ Therefore acceleration or deceleration are both possible. Physically the 1/γ root is related to the sub-linear dispersion ⟨r2 ⟩ ∼ tγ , whereas the provenance of f (p; γ) comes from the essence of the memory operator, as is seen from (2.16c). The former is a macroscopic spatial aftermath of the sluggishness induced by the sub-diffusive transport. The latter is quite unintuitive: seemingly long past evolution history asserts its influence via the convolution nature of (2.1e), or upon reduction by convolving u′ with a neither even nor odd function Dγy u in the second integral in (2.16c). At the limit of stagnation γ Ð→ 0+ said convolution outweighs the slow dispersion sufficiently to beget unbounded growth, voiding the validity of (2.16). In the original system (2.1) this is evinced by the pole of the Gamma function in (2.1e). Henceforth the attention is centred on the recovery scenarios, as they are responsible for a complete gamut of dynamics pertaining to system (2.16). 4.2. Instantaneous events: recovery sequences A limited assortment of solutions for sequences of the recovery type are reported to establish the main capabilities of system (2.16) as regards diverse drift behaviours subject to transient diffusion mechanisms. 4.2.1. Single spike Albeit the case of a single spike is degenerate from the aspect of asymmetry, it is instructive as regards basic expectations from a solution of (2.16) subject to the foregoing piecewise constant functions γ and D. In general the equilibrium constellation of n spikes is that of equidistantly distributed loci and identical height (with half the distance to the boundaries) 1/(βm−s−1) √ ⎛2 D 1 ⎞ 2i + 1 H (0) ∣ = tanh √ , i = {0, 1, . . . , n − 1}, (4.1) , xi ∣ = −1 + eq eq n ⎝ bm ⎠ n D whereby for a single spike x0 ∣ eq = 0. Figure 7 depicts the height and locus of the spike subject to the four recovery scenarios for γ and D (equation (3.2) and figure 5). Note the twofold discontinuity: the drift velocity x′i is discontinuous at moments (0) of cross-over in γ due to f (p; γ) by equation (2.16a), implying piecewise differentiable xi , whilst Hi is discontinuous at moments of cross-over in D by (2.16b) and (2.16d). By attentive scrutiny of panels 7(a) and 7(b), and juxtaposition with the transition scenarios some intriguing occurrences are gleaned. These are listed below. ‡ Albeit a cross-over in D abruptly adjusts the spike height, the effect on the drift velocity (slope of x0 on two sides of the transition point) might be virtually indistinguishable (seen from first two transitions in scenario D ⥯ γ ↟, thick black curves). The apparent monotonicity of the height response to changes in D is a most frequent occurrence, notwithstanding not to be deemed a universal inference. ‡ The spike height might be virtually irresponsive to discontinuities in γ (seen from the second transition in scenarios γ ⥯ D ↟ and γ ⥯ D ↡, thin curves). 8 Y. Nec Transitions in anomalous spike patterns ‡ Anomalous diffusion might tangibly speed up the motion, a phenomenon the more pronounced farther from the fixed point and responsible for markedly non-linear spike evolution (seen from the first transition in scenarios γ ⥯ D ↟ and γ ⥯ D ↡, thin curves). It must be noted that for more complex constellations the equilibration might not be monotonic in height, albeit remaining monotonic in the direction of drift as long as a stable regime is maintained. ‡ Despite at times similar impact of certain changes in the anomaly index γ and diffusion coefficient D, the exact sequence and cross-over timing determine the evolution history with the various stages being generally non-interchangeable. This simple case is a harbinger of phenomena to be expected for constellations of several spikes. Of most import is the evidence the reduced system (2.16) is capable of capturing variegated equilibration scenarios. Since each concentration spike corresponds to a cluster of molecules going through a diffusion and reaction process, it is to be inferred that one, transitions in γ and D together suffice to create a plethora of equilibration behaviours, and two, scrutiny of a given drift (0) path in {t, xi , Hi } space does not allow for an unequivocal illative determination of the input scenario. The analysis of the constellations below is done with these remarks in mind. 4.2.2. Two spikes For a pair of spikes by (4.1) the equilibrium state is {x0 , x1 } eq = { − 12 , 12 } and Hi (0) = H (0) ∣ , i = {0, 1}, the equilibrium height at each evolution stage conforming to the current diffusion coefficient. Whilst H (0) ∣ eq eq itself might change between consecutive stages within one transition sequence and the final value for various scenarios might not equal, the ratio (0) (0) H1 /H0 always approaches unity. This feature is useful in adjudgement of the proximity of a given pattern state to its equilibrium, along with x1 /x0 tending to −1. Figure 8 depicts the evolution of a two spike asymmetric pattern with transitions based on a recovery type scenario (equation (3.2) and figure 5). Panel 8(a) visualises how close to the equilibrium point the system is at the end of the simulation. Comparison of scenarios D ⥯ γ ↟ and D ⥯ γ ↡ (thick black and grey curves) reveals that anomalous diffusion speeds up the equilibration far from the fixed point, whilst in the vicinity thereof the impact of a cross-over in γ is minor. Overall a pattern starting as significantly anomalous reaches equilibrium faster, however such acceleration is not always monotonic over the entire range of γ. A detailed analysis of this non-monotonicity was presented in Nec [2018]. Panel 8(b) shows the relative location x1 /x0 versus time, all curves starting at − 31 and approaching −1. An easily discernible fact is that a cross-over in γ yields a markedly non-linear drift. An unrelated observation is that by juxtaposition of scenarios D ⥯ γ ↡ and γ ⥯ D ↡ (thick and thin grey curves), there is no correspondence between a pattern’s proximity to the fixed point as regards the height and locus: the pattern, whose evolution is shown by the thin grey curve, has the loci closer to equilibrium than those pertaining to the pattern with the thick grey path, whereas for the heights the contrary holds. Panels 8(c) and 8(d) give the phase portraits for both spikes, accentuating the fact that two not immediately comparable scenarios, D ⥯ γ ↟ and γ ⥯ D ↟ (thick and thin black curves), can result in nearly coincidental evolution. The initial stage differs due to the anomaly index, however between the first and second transitions the two curves are very close. After the third transition the path is nearly identical, only the initially anomalous pattern manifests a phase delay. These observations are by far not rare, thereby substantiating the difficulties faced by researchers interpreting experimental findings in systems with hybrid regular and anomalous diffusive processes. The numerical solution of (2.16b) forthwith becomes convoluted when the dimensionality increases, often presenting difficulty for higher values of the combination βm − s. This parameter is related to the pattern’s stability, cf. Iron et al. [2001]. Since the aim of the current study is a qualitative investigation, no particular perseverance was applied to design specific schemes that will successfully track the evolution very close to the onset of instability. Results are omitted whensoever improvement of an initial guess or further refinement of convergence and precision criteria in all numerical procedures did not yield graphically indistinguishable solutions within a feasible time frame. Addressing these numerical challenges is a topic of future studies. 4.2.3. Three spikes For a triple of spikes by (4.1) the equilibrium state is {x0 , x1 , x2 } eq = {− 32 , 0, 23 } and Hi (0) = H (0) ∣ , i = {0, 1, 2}. As before eq the equilibrium height itself might differ between the various stages of one scenario or respective parts of comparable transition sequences, however the height ratio of any pair out of the three spikes always approaches unity. Thus the ratios x1 /x0 and x1 /x2 and corresponding height ratios, tending to zero and unity respectively, were used to measure the proximity to equilibrium. Figure 9 gives a typical evolution of a three spike asymmetric pattern with transitions based on a recovery type scenario (equation (3.2) and figure 5). Phenomena described for constellations of fewer spikes recur, and some new features appear. A triple of spikes is the lowest dimensionality that allows for non-monotonic height ratio within the same evolution stage, i.e. with no cross-over in γ or D. Panels 9(a) and 9(b) exemplify this phenomenon for t×1 < t < t×2 . Comparison of panels 9(c) and 9(d) makes for an interesting observation: for the initial state used (identical for all solutions depicted in this figure) the ratio of loci of central to right spikes follows nearly the same course of evolution for two pairs of scenarios: {D ⥯ γ ↟, γ ⥯ D ↟ } and {D ⥯ γ ↡, γ ⥯ D ↡ }. Within each pair the transition sequences are by and large disparate, yet the location and hence velocity of the central and right spikes are very similar. Again, this is not 9 Y. Nec Transitions in anomalous spike patterns a rare occurrence, implying that being presented with only the drift path, inferring the correct transition scenario would be impracticable. From a physics vantage point this conforms to diversity in the sense that a quantitatively very similar outcome can be obtained via different imbedded diffusion mechanisms. 4.2.4. Four spikes For a quadruple of spikes by (4.1) the equilibrium state is {x0 , x1 , x2 , x3 } = { − 43 , − 41 , 14 , 34 } and Hi (0) = H (0) ∣ , i = {0, 1, 2, 3}. As evidenced in figure 10, typical evolution paths for this constellation appear similar to patterns of fewer spikes. The height ratio of the left pair of spikes, panel 10(a), strongly resembles the reciprocal ratio of central to right spikes in a triple, panel 9(b). The height ratio of the right pair of spikes, panel 10(b), is reminiscent of the ratio for a pattern of two spikes, 8(a). Loci ratios of central to outer spikes, panels 10(c) and 10(d), in addition to illustrating a compelling case of path similarity for disparate diffusive mechanisms, have all the hallmarks of respective central to outer ratios for a triple, panels 9(c) and 9(d). The phase portraits of the left pair, panels 10(e) and 10(f), are somewhat more diverse than e.g. 8(c), nonetheless presenting little qualitative novelty, whilst the right pair’s portraits, panels 10(g) and 10(h), are closely alike with that of the right spike out of a pattern of two, panel 8(d). In light of the above integrating (2.16) for patterns of higher dimensionality is unlikely to yield qualitatively new results. eq eq 4.3. Approximation of smooth transitions Figure 11 depicts the transition with 200 equal steps for each of γ(τ ) and D(τ ), the overall number of transitions dictated by the union set of cross-over moments. The number of steps was chosen so that the smooth cross-over behind the gradual passage {D , γ } Ð→ {D , γ } is easily envisioned, however the underlying discretisation still evident in regions of higher curvature. This is apparent in the height response, where discontinuity points conform to jumps in D alone. The locus x0 (τ ) is by definition continuous, but not differentiable at cross-over points, here the dearth of smoothness indiscernible upon visual scrutiny. Albeit formally the response of velocity x′0 (τ ) is discontinuous whensoever γ(τ ) or D(τ ) are, the fineness of this discretisation suffices to infer the fully smooth response merely by observation or simple numerical techniques for signal noise filtering, not applied here for the clarity of transition effects. Note that overall a transition in both γ(τ ) and D(τ ) occurs with forms (3.3) being symmetric about the same moment τ× , however at each cross-over moment only one parameter is discontinuous, rendering this scenario a refined extension of the instantaneous transition γ ↡ D ↡. Therefore panel 11(b) can be qualitatively compared e.g. to the second transition in panel 6(a) (bear in mind that the algebraic-differential system (2.16) is autonomous and hence translationally invariable in time). Similarly panel 11(c) can be juxtaposed with the first transition in 6(b). For completeness figure 12 shows similar results for a constellation of two spikes. Patterns with more spikes proved equally amenable to this technique and are not shown. 5. Discussion A system of two fractional partial differential equations governing the dynamics of two diffusing and reacting components was reduced by means of matched asymptotic expansions to describe the evolution of a constellation of spikes in a quasiequilibrium regime. The novel feature investigated herein is the possibility to capture transitions in both the diffusion coefficient D and anomaly index γ. The asymptotic description is valid for any order of unity variation of D(τ ), and strictly piecewise constant function γ(τ ) to uphold the matched asymptotic expansions leading to system (2.16), however is readily extended to generic functions γ(τ ) that are amenable to an approximation as a sum of weighted Heaviside step functions. The fineness of discretisation required for a visibly smooth cross-over was found to be unencumbering in terms of computation time and well proportioned to the gradients of the function γ(τ ). This is an analytical reduction that effectively simplifies variable order fractional partial differential equations for slowly varying D and γ. System (2.16) should be regarded as phenomenological, i.e. manifesting dynamics that transcend the specifics of particular parameter choices, thereby reproducing generic behaviour observed in unrelated physical systems. From a broader point of vantage this study recreates qualitative transient phenomena documented in recent experimental works in systems involving diffusion mechanisms with memory. Certain molecular events incur an abrupt shift in the diffusion coefficient, whilst not being classified as reactions, thus not entailing a necessity to deem the new configuration a distinct species. Amongst these might fall interactions with a catalyst affecting the molecule’s affinity with the medium fluid or tridimensional arrangement; a binding of a ligand or dissociation from one, whose impact on the modelled reaction is neutral, but affects the diffusive properties, such as a molecular motor or other progress aiding or hindering compound; or any other interaction that preserves the modelled reactive functionality, but alters the diffusive ability. On the other hand, a shift in the anomaly index might relate to a change in the medium, such as crossing into an essentially more or less crowded environment, or one, where the effective spatial dimension, possibly fractional, is distinct. Instances of this ilk are ingress into an organelle (nucleus, mitochondrion etc.) or egress therefrom, passage from a tri-dimensional diffusion within the cell to lateral diffusion within a membrane and more. In the Gierer-Meinhardt model a shift in the diffusion coefficient is tantamount to an adjustment of height, i.e. concentration of the reactant within the spike, conforming to a molecular cluster in more general models. A cross-over of the anomaly index begets a sudden acceleration or deceleration of the spike locus. The height and velocity are coupled. These effects were found to be non-monotonic or having a non-intuitive dependence on some of the parameters. It was not a rare occurrence to create a similar drift pattern by notably disparate diffusion scenarios. On the other hand, a loosely chosen and extremely limited – bearing in mind the dimensionality of the parameter space – set of transition sequences sufficed to span all qualitatively distinct outcomes. As this model is undeniably a fundamental and generic system as 10 Y. Nec Transitions in anomalous spike patterns regards the attributes of both diffusion and reaction processes, and in view of the fact that the foregoing findings did not require a dedicated search within the input parameter space, it stands to reason said findings transcend the particulars of the system used for this analysis. Illatively this evidence substantiates the difficulties experienced in interpretation of measurements and subsequent distinction between different diffusion regimes. If devising any number of input sequences, where diffusion changes from anomalous to normal and furthermore is sped or slowed due to shifts in the diffusion coefficient, with outcomes bearing very strong resemblance requires so little effort in a system as simple as this, obviously in a complex natural system backward inference of the transitions in diffusion mechanisms based on indirect measurements that in themselves are an intricate art, is an exceptionally difficult task. Indeed until the cross-over moments and diffusion mechanisms can be positively identified independently of the dispersion measurement, e.g. by understanding the physics of the system in question, it will remain a virtually impossible task. Appendix A. Fractional chain rule at the limit ǫ Ð→ 0 To prove the existence of a spike solution with sub-diffusion one requires an asymptotic form of the fractional time derivative for the inner solution Ai (yi (τ ), τ ) as pertains to the first argument, given by (2.1e) as: ∂τγ Ai (yi (τ )) = − 1 Γ (−γ) ∫ ⎧ ⎫ ⎪ x − xi (τ − ζ) ⎪ x − xi (τ ) ⎪ ⎪ dζ ⎨Ai ( ) − Ai ( ) ⎬ γ(τ )+1 . ⎪ ⎪ ǫ ǫ ζ ⎪ ⎪ 0 ⎩ ⎭ τ (A.1) Propose a new integration variable ξ ≡ (xi (τ − ζ) − xi (τ ))/ǫ. It is desired to invert the relation asymptotically when ǫ ≪ 1. The functions xi (τ ) and Ai (yi (τ ), τ ) are both infinitely differentiable under the consistency condition that adjacent spikes overlap no further than the exponentially small tails. Albeit this holds true with regular diffusion as well, the ordinary chain rule requires only the existence of the first derivative. Since xi ∈ C ∞ , the expansion xi (τ − ζ) = xi (τ ) − dxi 1 d2 x i 2 ζ+ ζ − +... dτ 2 dτ 2 (A.2) exists and is convergent. Then ξ = ǫ−1 (− Away from fixed points, where 1 d2 x i 2 dxi ζ+ ζ − + . . .) . dτ 2 dτ 2 (A.3a) dxi = 0, the series is inverted to give dτ ζ = −( 1 d2 x i 2 dxi −1 ) (ǫ ξ − ζ + − . . .) , dτ 2 dτ 2 dxi ≠ 0. dτ (A.3b) The series was never truncated, therefore (A.3b) is exact and holds everywhere bar the fixed points. Recursive substitution into higher powers of ζ gives dxi −1 ζ ∼ ǫ (− (A.3c) ) ξ + O (ǫ2 ) . dτ Thence for ǫ ≪ 1 equation (A.1) becomes ∂τγ Ai (yi (τ )) ∼ −ǫ−γ dxi dxi γ−1 γ ∣ ∣ Dyi Ai (yi ) , dτ dτ (A.4a) ⎧ ⎫ ⎪ ⎛ dxi ⎞⎪ ⎪ ⎪ dξ ⎨Ai (yi ) − Ai yi + sgn ( ) ξ ⎬ γ+1 . (A.4b) ⎪ dτ ξ ⎝ ⎠⎪ ⎪ ⎪ 0 ⎩ ⎭ Result (A.4) is a fractional equivalent of the chain rule for an infinitely differentiable function Ai in the limit ǫ Ð→ 0. dxi 1 Dγyi Ai (yi ) ≡ sgn ( ) dτ Γ (−γ) ∫ ∞ Acknowledgements. The author is grateful to J. Tzou (University of British Columbia) for advice on implementation of numerical schemes. References Bakalis, E., Höfinger, S., Venturini, A., 2015. Crossover of two power laws in the anomalous diffusion of a two lipid membrane. J. Chem. Phys. 142, 215102. Banks, D. S., Fradin, C., 2005. Anomalous diffusion of proteins due to molecular crowding. Biophys. J. 89, 2960–2971. ben Avraham, D., Havlin, S., 2000. Diffusion and reactions in fractals and disordered systems. Cambridge University Press, Cambridge. Bijeljic, B., Rubin, S., Sher, H., Berkowitz, B., 2011. Non-fickian transport in porous media with bimodal structural heterogeneity. J. Contaminant Hydrology 120-121, 213–221. Cabal, G. G., Genovesio, A., Rodrigues-Navarro, S., Zimmer, C., Gadal, O., Lesne, A., Buc, H., Feuerbach-Fournier, F., Olivo-Marin, J. C., Hurt, E. C., Nehrbass, U., 2006. SAGA interacting factors confine sub-diffusion of transcribed genes to nuclear envelope. Nature 441, 770–773. Dentz, M., Cortis, A., Sher, H., Berkowitz, B., 2004. Time behavior of solute transport in heterogeneous media: transition from anomalous to normal transport. Advances in Water Resources 27, 155–173. 11 Y. Nec Transitions in anomalous spike patterns Gierer, A., Meinhardt, H., 1972. A theory of biological pattern formation. Kybernetik 12, 30–39. Hahn, K., Kärger, J., 1995. Propagator and mean-square displacement in single-file systems. J. Phys. A: Math. Gen. 28, 3061–3070. 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. Klages, R., Radons, G., Sokolov, I. M., 2008. Anomalous transport: foundations and applications. Wiley-VCH, Weinheim. Metzler, R., Klafter, J., 2000. The random walk’s guide to anomalous diffusion: a fractional dynamics approach. Phys. Rep. 339, 1–77. Nec, Y., 2017. Dynamics of spike solutions in Gierer-Meinhardt model with time dependent diffusivity. J. Math. Anal. Appl., DOI http://dx.doi.org/10.1016/j.jmaa.2017.08.027. Nec, Y., 2018. Spike solutions in Gierer-Meinhardt model with a piecewise constant anomaly exponent. Comm. Nonlinear Science Num. Sim. 54, 267–285. 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. Ramirez, L. E. S., Coimbra, C. F. M., 2009. On the selection and meaning of variable order operators for dynamic modeling. International J. Diff. Eq. 2010, 846107. Sanabria, H., Kubota, Y., Waxham, M. N., 2007. Multiple diffusion mechanisms due to nanostructuring in crowded environments. Biophys. J. 92(1), 313–322. Saxton, M. J., 1994. Anomalous diffusion due to obstacles: a Monte Carlo study. Biophys. J. 66, 394–491. Saxton, M. J., 2007. A biological interpretation of transient anomalous subdiffusion. i. qualitative model. Biophys. J. 92(4), 1178–1191. 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. 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 Gierer-Meinhardt 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. Nonlinear Science 3, 209–264. Ward, M. J., Wei, J., 2005. The slow dynamics of two-spike solutions for the Gray-Scott and Gierer-Meinhardt systems: competition and oscillatory instabilities. SIAM J. Appl. Dyn. Sys. 4(4), 904–953. 12 Y. Nec Transitions in anomalous spike patterns −0.47 H0 x0 (0) 0.27 − 12 0.23 0 15 t 30 0 15 30 t (a) Spike locus and height (0) max H0 ∣ ǫ>0 (0) − H0 ∣ ǫÐ→0 0.01 0.005 0 0 0.1 ǫ (b) Initial height deviation Figure 2. Upper panel: comparison of single spike locus (left) and height (right) ensuing by systems (2.1) and (2.16), solid and dashed (time scaled τ = ǫ2 t) curves respectively, for a smooth D(τ ) and γ ≡ 1. Dotted curves: identical to solid, but with compensation for the finiteness of ǫ (see text for detail). Parameters for equation (3.3a): D = 2, D = 1, εD = 2, t× = 15. Initial state x0 (0) = − 21 . Other parameters: ǫ = 0.05, τ0 = 0.01. Chemical exponents used (p, q, m, s) = (2, 1, 2, 0). Qualitatively similar results obtained for sets (2,1,3,0), (3,2,2,0) and (4,2,2,0). Discretisation of (3.3a) into 200 equal steps yielded a solution visually indistinguishable from the dashed curve shown. Lower panel: max H0 ∣ (0) − H0 ∣ (0) , uncompensated for deviation of finite ǫ solution as given by (2.1) from initial quasi-equilibrium height at the limit ǫ Ð→ 0 1 versus ǫ (solid line), compared to 10 ǫ (dotted line). ǫ>0 13 ǫÐ→0 Y. Nec Transitions in anomalous spike patterns −0.47 H0 x0 (0) 0.27 0.23 − 12 0 15 t 30 0 15 30 t Figure 3. Comparison of single spike locus (left) and height (right) ensuing by systems (2.1) and (2.16), solid (dotted with compensation for finite ǫ error) and dashed (time scaled τ = ǫ2 t) curves respectively, for an instantaneous jump D(t < t× ) = 2, D(t ⩾ t× ) = 1, t× = 15, and γ ≡ 1. Initial state x0 (0) = − 12 . Other parameters: ǫ = 0.05, τ0 = 0.01. Chemical exponents used (p, q, m, s) = (2, 1, 2, 0). Qualitatively similar results obtained for sets (2,1,3,0), (3,2,2,0) and (4,2,2,0). γ↡ D↡ γ↡ D↟ γ D γ γ D D D γ γ γ D 0 t×1 D 0 t×2 t×1 t×2 t t γ↟ D↡ γ↟ D↟ γ γ D γ γ D D D γ γ D 0 t×1 D 0 t×2 t t×1 t×2 t Figure 4. Transition scenarios of the permanent change type, equation (3.1): anomaly index (left ordinate, diamonds) and diffusion coefficient (right ordinate, asterisks). Linestyles (colour and thickness) are preserved in subsequent simulations for ease of reference. 14 Y. Nec Transitions in anomalous spike patterns D ⥯ γ↟ D ⥯ γ↡ γ D γ γ D D D γ γ γ D 0 t×1 t×2 D 0 t×3 t×1 t×2 t×3 t t γ ⥯ D↟ γ ⥯ D↡ γ γ D γ γ D D D γ γ D 0 t×1 t×2 D 0 t×3 t t×1 t×2 t×3 t Figure 5. Transition scenarios of the recovery type, equation (3.2): anomaly index (left ordinate, diamonds) and diffusion coefficient (right ordinate, asterisks). Linestyles (colour and thickness) are preserved in subsequent simulations for ease of reference. 15 Y. Nec Transitions in anomalous spike patterns H0 (0) 0.28 0.19 0 t×1 t×2 t (a) Spike height x0 −0.28 − 34 0 t×2 t×1 t (b) Spike locus Figure 6. Single spike drift subject to permanent change type scenario (consult equation (3.1) and figure 4, scenario notation and curve styles conform thereto). Scenarios: γ ↡ D ↡ thick black curve, γ ↡ D ↟ thick grey curve, γ ↟ D ↡ thin black curve, γ ↟ D ↟ thin grey curve. Cross-over moments marked by diamonds. Initial state x0 (0) = − 43 . Chemical exponents used (p, q, m, s) = (2, 1, 2, 0). Qualitatively similar results obtained for sets (2,1,3,0), (3,2,2,0) and (4,2,2,0). 16 Y. Nec Transitions in anomalous spike patterns H0 (0) 0.28 0.19 0 t×1 t×3 t×2 t (a) Spike height x0 −0.27 − 43 0 t×1 t×3 t×2 t (b) Spike locus Figure 7. Single spike drift subject to recovery type scenario (consult equation (3.2) and figure 5, scenario notation and curve styles conform thereto). Scenarios: D ⥯ γ ↟ thick black curve, D ⥯ γ ↡ thick grey curve, γ ⥯ D ↟ thin black curve, γ ⥯ D ↡ thin grey curve. Cross-over moments marked by diamonds. Initial state x0 (0) = − 43 . Chemical exponents used (p, q, m, s) = (2, 1, 2, 0). Qualitatively similar results obtained for sets (2,1,3,0), (3,2,2,0) and (4,2,2,0). 17 Y. Nec Transitions in anomalous spike patterns (0) H1 /H0 (0) 0.97 0.66 0 t×1 t×3 t×2 t (0) (0) (a) Height ratio H1 /H0 x1 /x0 − 13 −0.83 0 t×1 t×3 t×2 t (b) Locus ratio x1 /x0 18 Y. Nec Transitions in anomalous spike patterns H0 (0) 0.17 0.16 − 43 x0 −0.55 (0) (c) Left spike phase portrait H0 (x0 ) H1 (0) 0.16 0.11 0.46 1 4 x1 (0) (d) Right spike phase portrait H1 (x1 ) Figure 8. Drift of two spike pattern subject to recovery type scenario (consult equation (3.2) and figure 5, scenario notation and curve styles conform thereto). Scenarios: D ⥯ γ ↟ thick black curve, D ⥯ γ ↡ thick grey curve, γ ⥯ D ↟ thin black curve, γ ⥯ D ↡ thin grey curve. Cross-over moments marked by diamonds. Initial state x0 (0) = − 34 , x1 (0) = 14 . Chemical exponents used (p, q, m, s) = (2, 1, 2, 0). Qualitatively similar results obtained for sets (3,2,2,0) and (4,2,2,0). 19 Y. Nec Transitions in anomalous spike patterns (0) H1 /H0 (0) 0.97 0.91 0 t×1 t×3 t×2 t (0) (0) (a) Central to left spike height ratio H1 /H0 (0) H1 /H2 (0) 1.14 1.03 0 t×1 t×3 t×2 t (0) (0) (b) Central to right spike height ratio H1 /H2 20 Y. Nec Transitions in anomalous spike patterns x1 /x0 1 4 0.12 0 t×3 t×2 t×1 t (c) Central to left spike locus ratio x1 /x0 x1 /x2 −0.16 −1 0 t×1 t×3 t×2 t (d) Central to right spike locus ratio x1 /x2 Figure 9. Drift of three spike pattern subject to recovery type scenario (consult equation (3.2) and figure 5, scenario notation and curve styles conform thereto). Scenarios: D ⥯ γ ↟ thick black curve, D ⥯ γ ↡ thick grey curve, γ ⥯ D ↟ thin black curve, γ ⥯ D ↡ thin grey curve. Cross-over moments marked by diamonds. Initial state x0 (0) = − 54 , x1 (0) = − 15 , x2 (0) = 51 . Chemical exponents used (p, q, m, s) = (4, 2, 2, 0). Qualitatively similar results obtained for sets (2,1,2,0) and (3,2,2,0). 21 Y. Nec Transitions in anomalous spike patterns (0) H1 /H0 (0) 0.99 0.91 0 t×1 t×3 t×2 t (0) (0) (a) Left central to outer spike height ratio H1 /H0 (0) H2 /H3 (0) 1.26 1.01 0 t×1 t×3 t×2 t (0) (0) (b) Right central to outer spike height ratio H2 /H3 22 Y. Nec Transitions in anomalous spike patterns x1 /x0 5 7 0.36 0 t×1 t×3 t×2 t (c) Left central to outer spike locus ratio x1 /x0 x2 /x3 0.32 − 13 0 t×1 t×3 t×2 t (d) Right central to outer spike locus ratio x2 /x3 23 Y. Nec Transitions in anomalous spike patterns H0 (0) 0.091 0.084 − 87 x0 −0.79 (0) (e) Leftmost spike phase portrait H0 (x0 ) H1 (0) 0.084 0.08 − 85 x1 (0) (f) Central left spike phase portrait H1 (x1 ) 24 −0.28 Y. Nec Transitions in anomalous spike patterns H2 (0) 0.082 0.064 − 81 0.22 x2 (0) (g) Central right spike phase portrait H2 (x2 ) H3 (0) 0.081 0.05 0.69 3 8 x3 (0) (h) Rightmost spike phase portrait H3 (x3 ) Figure 10. Drift of four spike pattern subject to recovery type scenario (consult equation (3.2) and figure 5, scenario notation and curve styles conform thereto). Scenarios: D ⥯ γ ↟ thick black curve, D ⥯ γ ↡ thick grey curve, γ ⥯ D ↟ thin black curve, γ ⥯ D ↡ thin grey curve. Cross-over moments marked by diamonds. Initial state x0 (0) = − 87 , x1 (0) = − 58 , x2 (0) = − 81 , x3 (0) = 38 . Chemical exponents used (p, q, m, s) = (2, 1, 2, 0). Qualitatively similar results obtained for sets (3,2,2,0) and (4,2,2,0). 25 Y. Nec Transitions in anomalous spike patterns 1 γ D 2 1 4 1 0 0 30 t (a) Piecewise constant approximation of smooth γ(t) (left ordinate, diamonds) and D(t) (right ordinate, asterisks) H0 (0) 0.246 0.215 0 t (0) (b) Spike height H0 (t) 26 30 Y. Nec Transitions in anomalous spike patterns x0 −0.424 − 34 0 t 30 (c) Spike locus x0 (t) x′0 0.474 0.188 0 t 30 (d) Drift velocity x′0 (t) Figure 11. Single spike drift for a piecewise constant approximation of a smooth cross-over. Parameters for equations (3.3): D = 2, D = 1, εD = 2, γ = 1, γ = 41 , εγ = 5, t× = 15, ǫ = 0.1, discretisation with 200 step-like transitions. Initial state x0 (0) = − 34 . Chemical exponents used (p, q, m, s) = (2, 1, 2, 0). Qualitatively similar results obtained for sets (3,2,2,0) and (4,2,2,0). 27 Y. Nec Transitions in anomalous spike patterns (0) H1 /H0 (0) 0.866 0.718 0 t (0) 30 (0) (a) Height ratio H1 /H0 x1 /x0 − 13 −0.645 0 t (b) Locus ratio 28 30 Y. Nec Transitions in anomalous spike patterns H0 (0) 0.168 0.166 −0.623 − 34 x0 (0) (c) Left spike phase portrait H0 (x0 ) H1 (0) 0.144 0.121 0.401 1 4 x1 (0) (d) Right spike phase portrait H1 (x1 ) Figure 12. Drift of two spike pattern for a piecewise constant approximation of a smooth cross-over. Parameters for equations (3.3): D = 2, D = 1, εD = 2, γ = 1, γ = 14 , εγ = 5, t× = 15, ǫ = 0.1, discretisation with 200 step-like transitions. Initial state x0 (0) = − 34 , x1 (0) = 14 . Chemical exponents used (p, q, m, s) = (2, 1, 2, 0). Qualitatively similar results obtained for sets (3,2,2,0) and (4,2,2,0). 29