Spike solutions in Gierer-Meinhardt model
with a time dependent anomaly exponent

1

2

Yana Nec

Thompson Rivers University, 805 TRU Way, Kamloops, British Columbia, Canada,

3
4

email: cranberryana@gmail.com

5

Abstract

6

Experimental evidence of complex dispersion regimes in natural systems, where the growth of the
mean square displacement in time cannot be characterised by a single power, has been accruing
for the past two decades. In such processes the exponent γ(t) in hr2 i ∼ tγ(t) at times might be
approximated by a piecewise constant function, or it can be a continuous function. Variable order
differential equations are an emerging mathematical tool with a strong potential to model these
systems. However, variable order differential equations are not tractable by the classic differential
equations theory. This contribution illustrates how a classic method can be adapted to gain insight
into a system of this type. Herein a variable order Gierer-Meinhardt model is posed, a generic
reaction– diffusion system of a chemical origin. With a fixed order this system possesses a solution
in the form of a constellation of arbitrarily situated localised pulses, when the components’ diffusivity ratio is asymptotically small. The pattern was shown to exist subject to multiple step-like
transitions between normal diffusion and sub-diffusion, as well as between distinct sub-diffusive
regimes. The analytical approximation obtained permits qualitative analysis of the impact thereof.
Numerical solution for typical cross-over scenarios revealed such features as earlier equilibration
and non-monotonic excursions before attainment of equilibrium. The method is general and allows
for an approximate numerical solution with any reasonably behaved γ(t).
Keywords: fractional differential equations, matched asymptotic expansions, variable order differential equations, numerical estimates of memory integrals

1. Background

9
10
11
12
13
14
15
16
17
18
19
20
21
22

23

Beginning circa 1980 measurements of dispersion in various natural systems documented
behaviour that could not be captured faithfully by a linear growth of the mean square
displacement hr2 i ∼ t. In time such systems were termed anomalous. Spanning numerous
scientific fields, the underpinning diffusive processes can be catalogued by basic transport
properties, often including the Fickian diffusion as a special limit (Codling et al., 2008; Eliazar
and Klafter, 2011). The sui generis nature of each anomalous process notwithstanding, wide
sub-classes have been identified (Klages et al., 2008). The present contribution addresses a
type of anomaly named sub-diffusion due to the sub-linear dispersion hr2 i ∼ tγ with 0 < γ < 1.
One way to obtain this kind of transport is to generalise the integer derivative in the diffusion
equation to one of a fractional order (Metzler and Klafter, 2000). Mathematically an integer

Preprint submitted to Elsevier

7
8

June 12, 2017

24
25
26
27
28
29
30
31
32
33

2

34
35

order partial differential equation for the concentration of components u with a constant
diffusion coefficient D and source f (u)
ut = D∆u + f (u),

36

38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70

t > 0,

(1a)

t > 0.

(1b)

will become a fractional partial differential equation
uγt = D∆u + f (u),

37

r ∈ Ω,

r ∈ Ω,

Over a certain range of time the dispersion r = |r| in the domain Ω will exhibit sub-diffusion
with a mean square displacement that grows in time with the exponent γ. Within this
framework the index γ need not be constant, albeit traditionally it has been fixed.
Whilst the integer derivative operator is indisputably unique, multifarious fractional
derivatives exist. For a constant γ the most classic definition appears in Oldham and Spanier
(1974), whereas more specialised operators can be found for instance in Elliot (1993) and
Chen et al. (2010), as well as in studies of specific systems quoted below. Variable order concominants have drawn some attention as well, cf. Naber (2004) and Ramirez and Coimbra
(2009). Recent studies bespeak the potential purport of the latter operators as a modelling
tool for systems evincing a time dependent dispersion power hr2 i ∼ tγ(t) .
One group of applications manifests distinct scales γ at well defined periods of time, i.e.
a piecewise constant approximation of γ(t) is a practical approach. An example thereof is
a time limited process in a finite domain, as undergone by molecules in biological media
(Reynolds, 2005; Bakalis et al., 2015), often with an auxiliary signalling reaction triggering
an essential change in the spatial structure of the diffusing molecules or local medium (Cabal
et al., 2006; Torreno-Pina et al., 2016). McKinley et al. (2009) construct a flexible model in
the field of rheology, and Meerschaert et al. (2013) present a more formal framework.
Another group of applications employs a continuous exponent γ(t). Sun et al. (2009) give
examples of numerical solutions with a linear function in t (as well as exponents depending on
variables other than time). Saxton (1997) examines the distribution of diffusion coefficients
in single particle tracking and demonstrates how statistically the apparent mean square
displacement hr2 i can be highly non-linear. The same author discusses a smooth cross-over
between two regimes with constant γ and compares to non-transient diffusion, connecting
the results to experimental measurement techniques (Saxton, 2007). At times γ and the
conforming order of the diffusion equation are noted to depend on system variables such
as concentration (Chen et al., 2013) or temperature (Morgan and Spera, 2001), however
ultimately when the mean square displacement is measured, the exponent γ will be timedependent.
In light of the above, variable order fractional equations promise a significant advance in
questions of correct interpretation of experimental measurements and juxtaposition of various
particle tracking techniques. To date studies presented numerical solutions that could be
compared to classic cases or fitted to experimental results, but offered little analytical insight.
An analytical solution to a paradigm system, explicating the temporal evolution subject to
transitions between disparate diffusive regimes and possible new behaviour thereupon, will

3
be a valuable step. This study focusses on a phenomenological reaction – diffusion system
with an aim to relax the restriction of a constant order temporal derivative (mean square
displacement growth power), thereby permitting the description of cross-over from regular
diffusion to sub-diffusion and between sub-diffusive regimes.
2. Gierer-Meinhardt model

∂tγ a = ǫ2γ axx − a +

ap
hq

− 1 < x < 1,

am
τo ∂tγ h = hxx − h + ǫ−γ s
h

ax (±1, t) = hx (±1, t) = 0 ,

t > 0,

73
74

76
77
78
79
80
81
82
83
84
85
86

(2a)
87

− 1 < x < 1,

t > 0,

(2b)

h(x, 0) = h0 (x) ,

(2c)

88

a(x, 0) = a0 (x) ,

where a(x, t) and h(x, t) are the activator and inhibitor concentrations, respectively. Here
0 < ǫ ≪ 1, τo > 0, the quadruple of reaction exponents (p, q, m, s) satisfies
q > 0,

72

75

The reaction – diffusion system in question was first proposed by Gierer and Meinhardt
(1972) and describes the interaction of two chemical species: an activator component that
encourages the reaction process and an inhibitor, suppressing the rate thereof. When the
diffusivity ratio between the components is asymptotically small with the reaction terms
appropriately scaled, the system possesses a solution in the form of a set of spikes, i.e.
localised spots of high concentration. With a fixed index γ these solutions were obtained by
Iron et al. (2001) for regular diffusion, γ = 1, and by Nec and Ward (2012) for sub-diffusion,
0 < γ < 1. Here it is proposed to endow the problem with time variability of the anomaly
index γ that is also the order of the partial differential equations, but in a way that still
admits of a solution.
On a finite, one-dimensional domain with no flux boundaries the unified model reads

p > 1,

71

m > 0,

s ⩾ 0,

p−1
m
<
,
q
s+1

89
90

(2d)

and the index γ ranges 0 < γ ⩽ 1. When γ = 1, system (2) comprises two conventional
partial differential equations. When 0 < γ < 1, the derivative of order γ is defined (for a
univariate function) as
Z t
1
f (t) − f (t − ζ)
dγ
f (t) = −
dζ,
0 < γ < 1,
(2e)
γ
dt
Γ(−γ) 0
ζ γ+1
wherein Γ denotes the Gamma function. Further detail on this operator can be found in
Elliot (1993) and Nec and Ward (2012). Since the integration is with respect to ζ, it is
forthwith admissible to extend the order to be a time dependent function γ(t). Some forms
of memory operators include differentiation with respect to t or imbed the delay directly into
the order through the dependence γ(t − ζ), therefore requiring a more careful generalisation,

91
92
93

94
95
96
97
98

4

108

cf. Ramirez and Coimbra (2009). Henceforth the anomaly index is deemed to be a function
of time, designated thus explicitly γ(t) or merely γ for simplicity.
For a constant γ the spike solution is constructed by the method of matched asymptotic
expansions (Iron et al., 2001; Nec and Ward, 2012). On an infinite domain a similar construction is possible (Nec, 2016b). One of the tacit assumptions of this classic method is that the
partial differential equations are of a fixed order, as are the concomitant asymptotic scales,
together admitting a reduction to a system of non-linear algebraic equations for the spikes’
heights and ordinary differential equations for the drift of loci. As shown hereinafter however, this reduction can be extended to the case of multiple transitions between a sequence
of arbitrary orders within the range 0 < γ ⩽ 1 at equally arbitrary cross-over moments.

109

3. A pattern of n spikes

99
100
101
102
103
104
105
106
107

110
111
112

Consider an n-tuple of spikes centred at a set of arbitrary loci xi , {i = 0, . . . , n − 1}.
Classically the outer solution for the activator away from the loci xi is the trivial quiescent
state a(x, t) ≡ 0, whereas for the inhibitor h(x, t) expand
h ∼ h(0) (x, t) + o(1).

(3)

.
Within a narrow region about the locus the inner spatial variable yi = (x − xi ) ǫγ(t) is
required to make the differential term on the right hand side in (2a) of order O(1). The
corresponding slow time variable τ = ǫα(t) t. The temporal scale α is determined below. The
inner asymptotic solutions are set as


(0)
(1)
γ(t)
−α(t)
Ai (yi , τ ) = a xi + ǫ yi , ǫ
τ ∼ Ai (yi , τ ) + ǫγ(τ ) Ai (yi , τ ) + · · ·
(4a)
def

113
114
115
116

117

118
119
120
121



(1)
(0)
Hi (yi , τ ) = h xi + ǫγ(t) yi , ǫ−α(t) τ ∼ Hi (yi , τ ) + ǫγ(τ ) Hi (yi , τ ) + · · · .

(4b)

The unusual feature about (4) is the time dependence of the scales γ and α, as conventionally
these are constant. It is imperative to substantiate when an expansion of this ilk is valid.
Obviously, it is valid when γ and α are constant. Observe that it is also valid for any step-like
function:
X n



o
γ(t) =
γj θ t − t× j−1 − θ t − t× j ,
j = 1, 2, . . . ,
(5)
j

122
123
124
125
126
127
128
129

wherein 0 < γj ⩽ 1 are constant anomaly indices, θ(t) denotes the Heaviside step function,
t× j are the corresponding cross-over times satisfying t× j > t× j−1 and t× 0 = 0 without loss of
generality. With γ(t) ≡ 1 the diffusion mechanism is the regular Fickian diffusion entailing
a linear growth of the mean square displacement in time hr2 i ∼ t, recovering the classic
model (Gierer and Meinhardt, 1972). With 0 < γ(t) < 1 the dispersion is sub-diffusive,
i.e. hr2 i ∼ tγ(t) . A plain, fixed order process is obtained with all γj equal. A regular to
anomalous transition is given by γj = 1, 0 < γj+1 < 1 and vice versa, whilst 0 < γj < 1,
γj 6= γj+1 yield transitions between distinct sub-diffusive regimes. At any cross-over moment

5
t× j the scales of the entire asymptotic series (4) will be instantaneously adjusted.
Albeit the non-linearity introduced into series (4) through the dependence ǫγ(t) forthwith
precludes an extension to functions of a more generic nature, it is hereby conjectured that the
spike solution will persist (albeit not be described by the series above) for smooth functions
such that the changes in γ(t) occur on a time scale much faster than τ . For instance for an
inner layer transition
γ(t) =

t − t× j γj+1 + γj
γj+1 − γj
tanh
+
2
ε
2

n−1
X

(0) βm−s

Hj

131
132
133
134
135

(6)

this requirement is mathematically expressed as ε ≪ ǫ−α , posing virtually no restriction
on the inner layer width ε since ǫ ≪ 1 and α > 0, whilst an inner layer exists as such
for ε ∼ O(1) or ε ∼ o(1). Thence illatively follows that the functional form of γ(t) is
encumbered but little, in particular any temporal variation on O(1) scale will not interfere
with the formation of the spike and its slow drift.
For the solution construction to hold formally it must be assumed that γ(t) is a piecewise constant function as in (5). However, any reasonably behaved function γ(t) might be
discretised into a set of step-like segments. The combination of solutions will remain valid
as long as t ∼ O(1) or higher, i.e. no attempt is made to seek a truly smooth evolution, or
to express that more technically, the orders of magnitude are preserved. Small steps in γ are
permissible as long as the restriction on t is observed. Figure 1 shows a possible discretisation
to approximate the smooth function (6).
For any given segment the classic reduction yields the following algebraic-differential
(0)
system (consult Nec and Ward (2012) for technical detail). The spikes’ heights Hi (leading
order inner solution for the inhibitor) are governed by a set of non-linear equations
(0)
Hi (τ ) = −bm

130

G(xi ; xj ) ,

136
137
138
139
140
141
142
143
144
145
146
147
148
149
150

(7a)

j=0

where G(x, xi ) is the Green’s function

 − csch 2 cosh(xi − 1) cosh(x + 1)
G(x; xi ) =

− csch 2 cosh(xi + 1) cosh(x − 1)

and

bm =

Z

151

x < xi
(7b)
x > xi
152

um dy ,

(7c)

R

is a factor computed from the homoclinic u(y)
u(y) =

p+1
2

!1/(p−1)

sech

153

2/(p−1)

p−1
y
2

!

(7d)

6

γ(t)

γj

γj+1
t× j

t× j+1
t

Figure 1: A valid discretisation example of a smooth cross-over.

154
155

156
157

that controls the spike’s shape. In particular the leading order inner solution for the activator
is given by
q
(0) β
(0)
β=
,
(7e)
Ai = Hi u(yi ),
p−1
whence Ai (yi , τ ) decays exponentially as |yi | −→ ∞ and is localised in the vicinity of the
spike locus xi . The leading order outer solution for the inhibitor reads
(0)

h (x, t) = −bm

n−1
X

(0) βm−s

Hi

G(x; xi ) .

(7f)

i=0

158
159

The spikes’ loci xi (τ ) drift on the slow time scale τ = ǫα t, α = γ + 1, according to the
differential system






n−1
γ−1

 X
βm−s
qbm f (p; γ) 1 (0) βm−s 
dxi dxi
(0)
−
+
G
(x
;
x
)
,
=
G
(x
;
x
)
+
G
(x
;
x
)
+
H
H
x i
j
x i
i
x i
i
j
(0)

dτ dτ
2 i
(p + 1)Hi 


j=0


j6=i

(7g)

7

160



f p; γ(t) ≡

Z

1
dxi
Dγy u(y) ≡
sgn
Γ(−γ)
dτ

up+1 dy

R

Z ∞(
0

,Z

u′ (y)Dγy u dy ,

(7h)

R

dxi
u(y) − u y + ξ sgn
dτ

161

!)

dξ
.
ξ γ+1

(7i)

With regular diffusion f (p; 1) = 2(p + 1)/(p − 1), equation (7h) being a fractional generalisation thereof. Operator (7i) ostensibly distinguishes between leftward and rightward drift,
thus breaking the symmetry of motion on either side of an equilibrium point. This is an
explicit aftermath of the memory inherent to all temporal fractional derivatives. Nonetheless
no actual asymmetry ensues, since the operator’s only appearance is in the second integral
in (7h), acting on the
with the odd function u′ (y), so that

 even function u(y), in conjunction
the factor f p; γ(t) does not depend on sgn x′i (τ ). The function Dγy u is neither even nor
odd for 0 < γ < 1 (consult figure 2). The functional dependence of f (p; γ) is given in figure
3. At a cross-over time t× j the value will be instantaneously adjusted by a jump between
respective values γj and γj+1 on the curve corresponding to the same value of the kinetic
exponent p. Bar the small interval of non-monotonicity for p = 2 (uppermost curve in figure
3), diminution in γ will always entail an increase in f (p; γ).
Albeit classically the foregoing asymptotic reduction (7) was devised for a constant γ, a
painstaking checking verifies its validity for a piecewise constant function γ(t). This result
is not altogether trivial, because whilst the inner spatial scale involves the compound ǫγ , the
slow time scale involves ǫα with α = γ + 1, whence a jump in γ cannot be forthwith absorbed
in ǫ. By and large system (7) permits a qualitative understanding of a variable order dynamical system of a fractional order γ(t): within the limitations of the asymptotic theory unifying
an arbitrary sequence of piecewise constant segments of γ(t), it shows analytically that the
variable order will not have an overwhelming impact on the nature of the solutions. Given
that from a mathematical point of vantage relaxing the constraint of a fixed order is devastating to the application of traditional solution and analysis tools, and at the outset might
lead to unpredictable results, this insight supports the physical soundness of variable order
models. Past computational works showed examples of this inference, e.g. Sun et al. (2009);
Chen et al. (2013), however with the Gierer-Meinhardt model a more formal substantiation
is possible. Moreover, owing to the generic nature of the piecewise constant discretisation,
it is safe to surmise that for most systems possessing asymptotic solutions of any kind the
same procedure will conform to a valid description for an arbitrary piecewise constant sequence of orders. In addition it is illustrated hereinafter for specific spike constellations that
a discretisation within the delineated limits of the presented theory yields a phase plane
trajectory remarkably smooth visually, i.e. a piecewise constant sequence provides a good
approximation without exacting the toll of a full numerical simulation. This last feature
pertains solely to the Gierer-Meinhardt model, of course, and cannot be extended illatively
to other dynamical systems.
System (7) is expected to be accurate by comparison to the full system (2) up to the

162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196

8

1.5

Dγy u

x′i <0

0

−1
−y∞

y∞

0
y

Figure 2: Function Dγy u with the fractional operator Dγy defined in (7i) for equally spaced 0.1 ⩽ γ ⩽ 0.9 (solid curves, neither
dxi
even nor odd functions) generalising u′ (dashed curve, odd function). The computation was performed for
< 0 and
dτ
p = 2 upon regularisation of (2e) via integration by parts twice. The infinite bound of the integral was truncated at 3y∞ .
Qualitatively similar results were obtained for 2 ⩽ p ⩽ 5. Cf. figure 1 of Nec and Ward (2012), computed to higher accuracy
here for expositional purposes.

197
198
199
200
201
202
203
204
205
206
207
208
209
210



order of the correction term in series (4), i.e. O ǫγ(τ ) . The concomitant error with regular
diffusion is O(ǫ). Similar reduced reaction– diffusion systems with regular diffusion have
been shown to manifest excellent agreement with the full partial differential system, insofar
as to predict correctly the onset of instability, cf. Sun et al. (2005) and Tzou et al. (2011).
In stark contrast to integer order equations, partial differential equations of a fractional
order, even without singular perturbations, require custom devised schemes, since conventional schemes have proved to baffle the regular notions as regards consistence, stability and
accuracy (Meerschaert and Tadjeran, 2004). These difficulties are caused by the presence of
memory, having at times unpredictable interference with other terms in the equations to be
solved. In addition, time fractional operators contain integrals with singular kernels. Without exception successfully implemented schemes circumvented this fundamental quandary
by regularising the operator. The modification involved either a change in the kernel, such
as Caputo’s derivative, several examples of which appear in Ramirez and Coimbra (2009),
or addition of an initial term (Henry and Wearne, 2000) or both (Sun et al., 2009). Defi-

9

f (p; γ)

9

6

4

3

0

1
γ

9
1
⩽ γ ⩽ 10
. Cross marks show the normal values
Figure 3: Factor f (p; γ) for equally spaced 2 ⩽ p ⩽ 5 computed for 10
+
2(p + 1)/(p − 1). f (p; γ) −→ ∞ at the stagnation limit γ −→ 0 . Cf. figure 3 of Nec and Ward (2012), computed to higher
accuracy here for expositional purposes.

nition (2e) is not regularised, its particular form essential to preserve basic equilibrium and
spectrum related properties of (2) as a dynamical system (more detail can be found in Nec
(2016b)). Thus a numerical solution of (2) is impossible with standard computational means.
Nonetheless, the error bound above can be used to calculate the maximal value of ǫ that
would ensure the solution given by (7) not exceed the corresponding error in the normal
system by taking min γ(τ ).
4. Quasi-equilibrium solutions

212
213
214
215
216

217

System (7) is the reduction of (2) at the limit ǫ −→ 0. Different constellations of spikes were
investigated to elucidate the effect of anomaly index transition on the drift toward equilibrium. Equation (7g) is a somewhat unconventionally written ordinary differential equation
that can be cast into a better amenable to numerical integration form by x′i (τ ) x′i (τ )

211

γ−1

=

218
219
220

221

γ

sgn x′i (τ ) x′i (τ ) , and thence x′i (τ ) can be computed by evaluating the right hand side of
(7g), whose sign determines sgn x′i (τ ). Generally at every integration step the non-linear
system (7a) must be solved, albeit some simple patterns permit an explicit solution. In

222
223
224

10

225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243

244
245
246

247

the classic Runge-Kutta procedure employed here or other predictor-corrector methods (7a)
must be solved for all intermediate predictor steps as well.
A more fundamental peculiarity of (7g) is that the integration variable τ relates to the
physical time t through τ = tǫγ(t)+1 , therefore changing scale at every cross-over point t× j ,
prescribed by equation (5). In praxis without lost of generality every
 characterised
 interval
(0)
– to the values
by a fixed value γj was integrated with τ initialised to zero and xi , Hi
reached at the end of the preceding interval. The full variable τ was made monotonic upon
completion of integration by addition of τ× j = t× j ǫγj +1 .
Infra solutions are presented for two groups of transition scenarios. The first category
examines a sequence of significant jumps in γ with the purpose to exemplify the extent of
abruptness in the pattern’s response with respect to relative velocity of the spikes and their
heights; demonstrate the variability of possible response behaviours; present non-intuitive
or non-monotonic response to certain parameters; identify typical mechanisms responsible
for pattern evolution that might be considered advantageous for a natural system, such as
a capability of fast equilibration or quick non-linear excursions before attainment of equilibrium. The second category pertains to the smooth cross-over. Function (6) was used with
the corresponding piecewise constant sequence in figure 1. Despite the coarseness of this
discretisation all evolution trajectories appeared absolutely smooth. Thus only a concise
sample is shown depicting this observation.
4.1. Single spike
With n = 1 equation (7a) is forthwith soluble for the spike height as a function of the locus
location x0

−1/(βm−s−1)
(0)
H0 = − bm G(x0 , x0 )
,
(8a)

whereby the differential equation (7g) can be integrated in closed form for γ = 1 to yield



cosh(2x
)
−
1
cosh(2x
(0))
+
1
0
0
−4qτ
sinh(2x0 )
1


= ln
+ cosh(2) ln 
(8b)
p−1
sinh(2x0 (0))
2
cosh(2x ) + 1 cosh(2x (0)) − 1
0

248
249
250
251
252
253

0

or numerically for 0 < γ < 1. Moreover, for any γ the sole fixed point x0 = 0 is stable,
allowing to replace sgn x′i (τ ) = − sgn xi . From (7g) it is readily inferred that the only effect
of γ is the speed wherewith the spike approaches equilibrium. What is less immediate is
whether a change in γ incurs magnification or diminution of said speed. The magnitude of
the right hand side of (7g) is below unity for any practical values of parameters, exclude
for instance extremely small γ since lim+ f (p; γ) = ∞. Thereby x′i (τ ) will involve this
γ−→0

254
255
256
257
258

quantity’s 1/γ > 1 power for any 0 < γ < 1. Bar a narrow region just below γ = 1 for p = 2,
f (p; γ) is a descending function of γ (figure 3). Therefore the 1/γ > 1 power and growth of
f (p; γ) are opposing influences, begetting the following behaviour: for a small increase in f ,
i.e. γ slightly below unity, the spike will slow down somewhat, for more significant growth of
f it will speed up, whilst for small γ its progress will again be delayed at the cross-over point

11
t× , this latter evolution occurring over a very narrow interval of γ. This is best exemplified
when f depends monotonically on γ, i.e. p > 2, although it is always correct. Figure 4 shows
a typical example of this phenomenon.

259
260

x′0 (τ )

max x′0

0

0

τ×
τ

Figure 4: Anomaly dependent acceleration and delay phenomenon during the drift of a single spike: velocity x′0 (τ ) for a
transition scenario γ(t < t× ) = γ1 , γ(t > t× ) = γ2 with γ1 = 1 and γ2 = 1 (dotted), γ2 = 12 (dash-dotted), γ2 = 0.0995
(dashed) and γ2 = 0.099 (solid). Note the non-monotonicity due to growth of f (p; γ) with diminution of γ counteracted by
1
the 1/γ power in equation (7g). Other parameters used: (p, q, m, s) = (3, 2, 2, 0), t× = 10, ǫ = 10
, x0 (0) = − 21 . Qualitatively
similar results obtained with (p, q, m, s) = (2, 1, 2, 0), (2, 1, 3, 0), (4, 2, 2, 0).
261

With a single spike, or more generally any constellation with only one drift equation
(more examples to be discussed hereinafter), the anomaly dependent
acceleration
or delay


(0)
th
plane. The slow
do not affect the trajectory traversed by the i spike in the xi , Hi
time variable τ parametrises this trajectory, hence transitions between different values of γ
impact solely the velocity of the spike motion, as visualised in figure 5. In light of the above
a general conclusion follows that for this degenerate pattern additional transitions will not
result in worthy of investigation phenomena.
4.2. Pair of spikes
A well explored constellation is the symmetric pattern, where the premise of reflective sym(0)
(0)
metry x0 < 0, x1 (τ ) ≡ −x0 (τ ) entails equal heights H0 (τ ) ≡ H1 (τ ) (Iron et al., 2001;
Nec and Ward, 2012; Nec, 2016a). System (7) effectively reduces to one forthwith soluble

262
263
264
265
266
267
268

269
270
271
272

12

298

non-linear equation and one ordinary differential equation, whose closed form solution for
γ = 1 is very similar to (8) and omitted here. Illatively from the single spike case, this
pattern is not interesting as regards transient phenomena due to variable anomaly indices,
as the trajectories are invariant albeit the parameterisation differs. Therefore it behoves one
to centre the attention on asymmetric arrangements, whereat system (7) comprises two pairs
of coupled equations.
The effect of varying anomaly index is better evident when initially the loci are positioned
significantly away from the symmetric constellation. Figure 6 shows the trajectories of
(0)
Hi (xi ), i = {0, 1} for a scenario of five distinct anomaly indices. Note the non-monotonicity
of the left spike’s height adjustment: with diminishing γ in the first two transitions the spike
grows higher, then lower as γ decreases further at the third transition, followed by slight
growth with γ returning to normal. The response of the right spike does not manifest obvious
conforming evolution, as it is located farther from its equilibrium point, its motion centred on
attainment thereof. In particular the left spike overshoots the equilibrium height for a spell to
compensate for the right spike’s position. Indeed eventually the two spikes are similarly close
to the respective equilibria (compare locations at trajectory ends), albeit never symmetric.
Juxtaposition with the propagation without anomaly reveals a twofold impact: the left spike,
initially situated near its equilibrium point, would have achieved the equilibrium height very
slowly but for the transition, whilst the right spike would have been nowhere close to the
equilibrium position or height within the same time frame. Thus the transition to anomaly
might be a powerful mechanism to speed up the pattern’s equilibration.
With a set of kinetic exponents corresponding to higher values of β, and thence a stronger
non-linearity, it is also possible to obtain a trajectory that includes an excursion away from
equilibrium before converging thereto. This occurs when the initial heights of the two spikes
are on different sides of the equilibrium value. Nonetheless here too the equilibration is
essentially faster with anomaly. Figure 7 depicts a situation of this ilk.

299

4.3. Triple of spikes

273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297

300
301
302
303
304
305
306
307
308
309
310
311
312


 

The symmetric arrangement x0 (τ ), x1 (τ ), x2 (τ ) ≡ x0 (τ ), 0, −x0 (τ ) again effectively reduces system (7)
equation of the autonomous type, to wit no differences

 to a single differential
(0)
in trajectories xi (τ ), Hi (τ ) will ensue regardless of the parameterisation in τ . Thus, as
before, only asymmetric constellations are of interest. Figure 8 gives a typical equilibration
trajectory. Interesting phenomena to observe are earlier approach towards the equilibrium
(0)
height Heq for all three spikes with a transition to sub-diffusion, possible overshooting of
(0)
Heq whilst drift with regular diffusion begets no such occurrence (right spike x2 ), and fur(0)
ther excursion from Heq by comparison to the regular trajectory (central spike x1 ). Figure
9 depicts the evolution for a combination of kinetic parameters that incur initial height on
(0)
different sides of the equilibrium value Heq . With this set of parameters, corresponding to
a higher value of β, it is possible to obtain an overshooting of equilibrium by the normal
trajectory, but not the anomalous one (left spike, x0 ), as well as a considerably delayed
overshooting (right spike, x2 ).

13
4.4. Quadruple of spikes
For any n > 3 a symmetric arrangement of either even


x0 (τ ), x1 (τ ), . . . xn/2−1 (τ ), −xn/2−1 (τ ) . . . , −x1 − x0

313

or odd



x0 (τ ), x1 (τ ), . . . x(n−1)/2−1 (τ ), 0, −x(n−1)/2−1 (τ ) . . . , −x1 − x0
 
number of spikes system (7) will require the solution of n2 non-linear equations for the
spike heights and an equal number of respective drift differential equations.
In this sense


(0)
a pattern of four spikes is the first case, where a non-smooth trajectory xi (τ ), Hi (τ ) is
expected with a symmetric initial state. Figure 10 shows such an example.
With an asymmetric initial state similar phenomena are obtained as with n ⩽ 3: considerable excursions, faster occurring changes in spike height and marked non-monotonicity.
Rather consistently the spikes drift a larger distance than subject to an uninterrupted normal diffusion regime. A typical example is given in figure 11. In light of these observations
constellations of n > 4 are expected to introduce more variability within the available degrees of freedom, but reveal no qualitatively new phenomena. It is noteworthy that the
same set of kinetic exponents that for pairs and triples engendered less typical evolution,
with n = 4 failed to yield an equilibrated pattern (with initial conditions and all transition
parameters identical to other sets used). This occurrence was further investigated and is
discussed hereinafter.
Figure 12 illustrates the smoothness of phase plane trajectories for the same constellation
and initial conditions, but a discretised transition sequence approximating a smooth crossover in figure 1. From the result it is seen that a full solution would not produce tangibly
smoother curves, proving the soundness of the approximation procedure. The computation
was performed for all foregoing constellations with a similar aftermath. The particular
example shown was chosen due to well defined non-monotonicity in height evolution that
was nonetheless traced smoothly.
5. Existence of solutions

314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334

335

(0)

System (7) is valid as long as Hi > 0 ∀i and all xi are sufficiently far apart. Being a
reduction of the full partial differential system (2) on the slow time scale τ , (7) cannot track
behaviour that breaches the assumptions of quasi-equilibrium.
One such event of relevance here is a termination of a solution branch due to a bifurcation.
To the best of the author’s knowledge this occurrence has not been formally explored in the
algebraic-differential system associated with the spike solution in Gierer-Meinhardt model,
albeit it has been observed in numerical
(Tzou, 2016). A constellation following

 simulations
(0)
a branch in the 2n-dimensional space xi , Hi , i = {0, . . . , n − 1} that terminates before
the equilibrium point is reached, will not be able to equilibrate. Only simulation of the
full partial differential system (2) can elucidate the shift to another solution branch and

336
337
338
339
340
341
342
343
344
345

14

370

possible changes in the constellation, a study beyond the ambit of the current analysis.
The existence
and concomitant termination of the trajectories in each of

 of a bifurcation
(0)
is an inherent property of the reduced Gierer-Meinhardt system and
the planes xi , Hi
not the result of sub-diffusion or transitions between distinct values of the anomaly index,
however a bifurcation might induce an earlier termination in the presence of anomaly or
beget unbounded growth. Typical examples appear hereunder, and no attempt is made to
map the parameter space effecting a bifurcation.
The first example is for a system with regular diffusion and no transitions, i.e. γ ≡ 1.
(0)
Figure 13 shows the spike velocity x′i (τ ) and trajectories Hi (xi ) for a quadruple of spikes.
Whilst the velocity is finite, note that the slope at the terminus point τtrm becomes infinite, as
d (0)
does
H . Nevertheless x′i (τ ) remains finite, as this is not a singularity of (7g). One must
dτ i
bear in mind the terminus is indeed attained – the integration ceases because (7a) possesses
no solution, in contrast to the equilibration itself with its exponentially or algebraically slow
approach towards the fixed point, or a finite time blow up, where computation becomes
impractical.
If γ is moderately small, a constellation might reach a branch terminus, whilst with γ ≡ 1
no bifurcation is encountered. An example of this is given in figure 14. This constellation
is initialised identically to the one used in figure 11 with sets of kinetic exponents begetting
no solution existence issues.
An unimpeded growth of x′i (τ ) might ensue for small values of γ, when the right hand
γ
side of (7g) exceeds unity in magnitude. Re-write the left hand side of (7g) as sgn x′i x′i ,
whence it is readily inferred that 1/γ begets a power growth. For instance, with the same set
of parameters as in figure 14, if another cross-over moment is added anywhere in the range
t× 1 < t× 2 < ttrm and γ3 < γ2 , the solution grows beyond the computer’s handling ability
within a few integration steps after the transition.

371

6. Discussion

346
347
348
349
350
351
352
353
354
355

356
357
358
359
360
361
362
363
364
365
366
367
368
369

372
373
374
375
376
377
378
379
380
381

The infinitesimally narrow spike solution, existing for the Gierer-Meinhardt model with
both regular diffusion (γ = 1) and sub-diffusion (0 < γ < 1), persists when the index
γ, traditionally a constant, is a piecewise constant function. To wit, the validity of the
asymptotic spike pattern extends to the situation when the system undergoes transitions
between regimes, whose mean square displacement evolution hr2 (t)i ∼ tγ is characterised by
distinct exponents. The number of transition points as well as their cross-over timings t× j
and corresponding anomaly indices 0 < γj ⩽ 1 are arbitrary.
At each transition γj −→ γj+1 occurring at t× j , j = 1, 2, . . ., the drift velocity x′i (τ ) has a
(0)
jump discontinuity, therefore the spike height Hi (τ ) and locus position xi (τ )are continuous,

(0)
but not smooth. Generally the non-smoothness is retained in the plane xi , Hi , with
(0)

382
383

the exception of trajectories Hi (xi ) where τ is an arbitrary parameterisation variable, i.e.
for those constellations that are effectively described by a sole non-linear equation for the

15
height and respective differential equation for locus drift, namely a single spike, a pair of
symmetric spikes, and a triple with central spike at equilibrium x1 ≡ 0 and outer spikes
moving symmetrically. Thence for n ⩽ 3 only asymmetric initial states yield qualitatively
distinct trajectories for a sequence of transitions by comparison to a constant γ(t). For
n ⩾ 4 both symmetric and asymmetric initial states give non-smooth trajectories in the
locus position – spike height planes. Near the equilibrium points transitions occasion but
small differences in velocity, therefore often the non-smoothness is virtually indiscernible.
The manifold alterations in the trajectory form range from minor and quantitative to
qualitatively novel behaviour. Depending on the specific parameters chosen, a pattern might
equilibrate significantly faster than with a single fixed anomaly index, be it normal or subdiffusive; the spikes might drift a larger distance in the process; the spike height might reach a
local maximum or minimum that would not have been attained without a series of transitions.
Recollecting that the spike height conforms to reactant concentration in the original chemical
model, these results bear on recent experimentally documented concepts related to reactions
within such complex biological media as cell membranes, organelles and nuclei. The most
prominent aftermath of permitting the conventionally constant anomaly index to vary in
time is the instantaneous adjustment of the spike’s velocity. This behaviour qualitatively
corresponds to a number of scenarios. One is a concatenation of events, when a cascade
of consecutive periods of molecule diffusion and subsequent reactions trigger each other in
turn, each with its own physical characteristics due to the locality of each drifting cluster and
ligands. When out of a set of identical diffusing molecules a sufficient number is bound at the
target, the concentration of that reactant diminishes, signalling the next reaction in the chain
and possible slowing down of the initial diffusion. Another is that docking of the designated
molecules will trigger a fast removal of the remaining reactant, thereby requiring a significant
acceleration of their spatial propagation. Yet another is diffusion through distinct media, like
penetration in and out of a nucleus or a cell or within such systems as the lymphatic vessels
and nodes, where multifarious stages can be characterised by different anomaly indices.
The findings presented are based on an asymptotic reduction of a system of partial
differential equations of a fractional order. Ideally it is desirable to obtain a numerical
verification in those parts of the parameter space, where the asymptotic solution is valid.
Numerical solution of partial fractional differential equations is notoriously difficult and
often defies conventional concepts of numerical analysis as regards consistence, accuracy and
stability of schemes (Meerschaert and Tadjeran, 2004). Here the challenge due to the kernel
singularity of (2e) is further compounded by an increase in the spike width with diminishing
γ, whereby the interval of feasible values of ǫ is narrower than the concomitant expected for
normal diffusion and integer order equations. Development of a suitable numerical scheme
is an enduring problem and topic of future research. All past successful schemes resorted
to a regularisation of the fractional operator, a course of action impossible here, given that
the operator is needed in its current form for the spike pattern to exist as well as certain
spectral properties to be retained.
The Gierer-Meinhardt system is a paradigmatic model, representing no particular natural
application, but capable of reconstructing general observable phenomena. Whist having once

384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425

16

448

more proved a powerful pattern forming tool, the limitation of the solution obtained here is in
the requirement that the anomaly index remain constant between any two transition points.
Whilst a formally continuous variation of γ might be a curious topic for future research,
the posed theory was shown to yield evolution trajectories appearing visibly smooth when a
smooth function γ(t) was discretised into a sequence of step-like segments. The limitations
on the validity of the presented asymptotic solution are not burdensome, therefore such
discretisation is feasible for most reasonably behaved γ(t). Moreover, the method itself is
generic and might be easily adapted to unrelated dynamical systems, as long as an asymptotic
solution with fixed scales is available. The degree of smoothness in the solution will depend
on the particular system at hand. It is furthermore conjectured that for γ(t) of an abrupt,
yet continuous variation, for instance an inner layer as in (6) with ε ≪ 1, away from the
transition core the solution obtained herein will faithfully capture the dynamics, whilst in
the vicinity thereof some delay in adjustment of the spike velocity will ensue, depending
on the specific dependence of γ(t). Perhaps an asymptotic solution might be furnished for
certain forms of γ(t) based on an independent small parameter describing the abruptness of
transition. With an unrestricted variation, i.e. γ(t) such that no regions can be approximated
by a constant, a disparate method of solution will be required, as the current asymptotic
approach becomes unfeasible.
The Gierer-Meinhardt system is further amenable to a generalisation involving a fractional Laplacian (spatial fractional derivative), thereby admitting solutions of spikes that
propagate according to Lévy flights, a type of super-diffusion (Nec, 2012). Extension of the
current analysis to include transitions between Lévy flights with distinct anomaly indices as
well as between sub-diffusive to super-diffusive regimes is another future problem.

449

References

426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447

450
451

452
453
454
455

456
457

458
459

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.
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.
Chen, W., Sun, H., Zhang, X., Korošak, D., 2010. Anomalous diffusion modeling by fractal
and fractional derivatives. Comp. Math. Appl. 59, 1754–1758.
Chen, W., Zhang, J., Zhang, J., 2013. A variable-order time-fractional derivative model for
chloride ions sub-diffusion in concrete structures. Frac. Calc. Appl. Analysis 16, 76–92.

461

Codling, E. A., Plank, M. J., Benhamou, S., 2008. Random walk models in biology. J. R. Soc.
Interface 5, 813–834.

462

Eliazar, I., Klafter, J., 2011. Anomalous is ubiquitous. Annals of Physics 326, 2517–2531.

460

17
Elliot, D., 1993. An asymptotic analysis of two algorithms for certain Hadamard finite-part
integrals. IAM J. Num. Analysis 13, 445–462.
Gierer, A., Meinhardt, H., 1972. A theory of biological pattern formation. Kybernetik 12,
30–39.
Henry, B. I., Wearne, S. L., 2000. Fractional reaction–diffusion. Physica A 276, 448–455.
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.
McKinley, S. A., Yao, L., Forest, M. G., 2009. Transient anomalous diffusion of tracer
particles in soft matter. J. Rheol. 53, 1487.
Meerschaert, M. M., Nane, E., Vellaisamy, P., 2013. Transient anomalous sub-diffusion on
bounded domains. Proc. Amer. Math. Soc. 141, 699–710.
Meerschaert, M. M., Tadjeran, C., 2004. Finite difference approximations for fractional advection – dispersion flow equations. J. Comput. Appl. Math. 172, 65–77.
Metzler, R., Klafter, J., 2000. The random walk’s guide to anomalous diffusion: a fractional
dynamics approach. Physics Reports 339, 1–77.
Morgan, N. A., Spera, F. J., 2001. A molecular dynamics study of the glass transition in
CaAl2 Si2 O8 : thermodynamics and tracer diffusion. American Mineralogist 86, 915–926.
Naber, M., 2004. Distributed order fractional sub-diffusion. Fractals 12, 23.
Nec, Y., 2012. Spike-type solutions to one dimensional Gierer-Meinhardt model with Lévy
flights. Studies in Appl. Math. 129, 272–299.
Nec, Y., 2016a. Dynamics of spike solutions in Gierer-Meinhardt model with time dependent
diffusivity. preprint.
Nec, Y., 2016b. Explicitly solvable eigenvalue problem and bifurcation delay in sub-diffusive
Gierer-Meinhardtmodel. 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.

463
464

465
466

467

468
469

470
471

472
473

474
475

476
477

478
479

480
481

482

483
484

485
486

487
488

489
490

Oldham, K. B., Spanier, J., 1974. Fractional Calculus. Academic Press, New York.

491

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.

493

492

18

494
495

496
497

498
499

500
501

502
503
504

Reynolds, A. M., 2005. On the anomalous diffusion characteristics of membrane-bound proteins. Physics Letters A 342, 439–443.
Saxton, M. J., 1997. Single-particle tracking: the distribution of diffusion coefficients. Biophys. J. 72, 1744–1753.
Saxton, M. J., 2007. A biological interpretaton of transient anomalous subdiffusion. i. qualitative model. Biophys. J. 92, 1178–1191.
Sun, H., Chen, W., Chen, Y. Q., 2009. Variable-order fractional differential operators in
anomalous diffusion modeling. Physica A 388, 4586–4592.
Sun, W., Ward, M. J., Russell, R., 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.

507

Torreno-Pina, J. A., Manzo, C., Garcia-Parajo, M. F., 2016. Uncovering homo- and heterointeractions on the cell membrane using single particle tracking approaches. J. Phys. D:
Appl. Phys. 49, 104002.

508

Tzou, J. C., 2016. private communication.

505
506

509
510
511

Tzou, J. C., Bayliss, A., Matkowsky, B. J., Volpert, V. A., 2011. Stationary and slowly
moving localised pulses in a singularly perturbed Brusselator model. Euro. J. Appl. Math.
22(5), 423–453.

(0)

x0 (τ )

H0 (τ )

19

(0)

H0 (τ× )
x0 (τ× )

(0)

H0 (0)

x0 (0)
τ×

τ×

0

τ

τ

H0

(0)

0

(0)

H0 (τ× )

(0)

H0 (0)
x0 (0)

x0 (τ× )

x0

Figure 5: Trajectory parameterisation example for a single spike and transition scenario γ(t < t× ) = 1, γ(t > t× ) = 15
(0)

(solid curves) and drift with no transition, γ ≡ 1 (dotted curves): locus position x0 (τ ) (upper left panel), H0 (τ ) (upper
(0)
right panel) and H0 (x0 ) (bottom panel).

Diamonds mark the cross-over moment τ× = t× ǫγ(t)+1 on all panels. Circle
marks the end of the indistinguishable dotted curve on the bottom panel: from the upper panels it is obvious that the spike
speeds up upon transition, therefore the solid trajectory (with transition) extends further than the dotted one (no transition).
1
, x0 (0) = − 43 . Qualitatively similar results obtained with
Other parameters used: (p, q, m, s) = (3, 2, 2, 0), t× = 50, ǫ = 10
(p, q, m, s) = (2, 1, 2, 0), (2, 1, 3, 0), (4, 2, 2, 0).

20

(0)

H0

(0)

Heq

(0)

H0 (0)
-0.42
x0

x0 (0)

(0)

H1

(0)

Heq

(0)

H1 (0)
x1 (0)

0.42
x1

Figure
of five intervals
characterised by

o parameterisation example for a pair of spikes and transition
 scenario
 
n 6: Trajectory
3
1
1
. Equilibrium position
, 10
, 1 with cross-over moments t× = {1, 2, 6, 8}. Initial position is x0 , x1 = − 52 , 10
γ = 1, 41 , 20




is x0 , x1 = − 21 , 12 . Equilibrium height, equal for both spikes, is marked by the dashed line. Dotted lines show the
trajectory for γ ≡ 1 (in upper panel continuing beyond the frame given). Diamonds mark the cross-over moments t× j . Other
1
parameters used: (p, q, m, s) = (3, 2, 2, 0), ǫ = 10
. Qualitatively similar results obtained with (p, q, m, s) = (2, 1, 2, 0), (4, 2, 2, 0).

21

(0)

H0

(0)

max H0

(0)

Heq

x0 (t× 4 )

x0

x0 (0)

(0)

H0

(0)

max H0

(0)

H0 (0)
x0 (t× 3 )

x0 (t× 2 )

x0

x0 (t× 1 )

22

(0)

H1

Heq

(0)

H1 (0)
x1 (t× 1 )

x1 (t× 4 )

x1

Figure
of non-monotonic equilibration for a pair of spikes with transition scenario of five intervals characterised by
n 7: Example o
3
1
, 10
, 1 with cross-over moments t× = {1, 2, 6, 8}: left spike exhibits trajectory excursion away from equilibrium
γ = 1, 41 , 20
prior to convergence (upper
enlarged in central 
panel), right spike’s 
response is monotonic (lower

 panel, height
 peak vicinity
panel). Initial position is

x0 , x1

=

1
. Equilibrium position is
− 52 , 10

x0 , x1

=

− 12 , 21 . Equilibrium height, equal

for both spikes, is marked by the dashed line (observe left spike begins above equilibrium height, whilst right spike is below).
Dotted lines show the trajectory for γ ≡ 1 (in upper and central panels continuing beyond the frame given). Diamonds mark
1
the cross-over moments t× j . Other parameters used: (p, q, m, s) = (2, 1, 3, 0), ǫ = 10
.

23

(0)

H0

(0)

Heq

(0)

H0 (t× 1 )
x0 (t× 4 )

x0

x0 (t× 1 )

(0)

H1

(0)

max H1

(0)

Heq

(0)

H1 (t× 1 )
x1 (t× 4 )

x1

x1 (t× 1 )

24

(0)

H2

(0)

max H2

(0)

Heq

(0)

H2 (t× 1 )
x2 (t× 1 )

x2

x2 (t× 4 )

Figure
by

o parameterisation example for a triple of spikes and transition scenario of
 five
 intervals characterised
n 8: Trajectory
3
1
1
3
. Equilibrium
, 10
, 1 with cross-over moments t× = {1, 2, 6, 25}. Initial position is x0 , x1 , x2 = − 35 , 10
, 10
γ = 1, 41 , 20


 
position is x0 , x1 , x2 = − 32 , 0, 23 . Equilibrium height, equal for all spikes, is marked by the dashed line. Dotted lines show
1
the trajectory for γ ≡ 1. Diamonds mark the cross-over moments t× j . Other parameters used: (p, q, m, s) = (3, 2, 2, 0), ǫ = 10
.
Qualitatively similar results obtained with (p, q, m, s) = (2, 1, 2, 0), (4, 2, 2, 0).

25

(0)

H0

(0)

Heq

(0)

H0 (t× 1 )
x0

x0 (t× 4 )

x0 (t× 1 )

(0)

H1

(0)

max H1

(0)

Heq

x1 (t× 4 )

x1

x1 (t× 1 )

26

(0)

H2

(0)

max H2

(0)

Heq

(0)

H2 (t× 1 )
x2 (t× 1 )

x2

x2 (t× 4 )

Figure
by
n 9: Trajectory

o parameterisation example for a triple of spikes and transition scenario of
 five
 intervals characterised
3
1
1
3
γ = 1, 41 , 20
. Equilibrium
, 10
, 1 with cross-over moments t× = {1, 2, 6, 25}. Initial position is x0 , x1 , x2 = − 35 , 10
, 10


 
position is x0 , x1 , x2 = − 32 , 0, 23 . Equilibrium height, equal for all spikes, is marked by the dashed line. Dotted lines show
1
the trajectory for γ ≡ 1. Diamonds mark the cross-over moments t× j . Other parameters used: (p, q, m, s) = (2, 1, 3, 0), ǫ = 10
.

27
(0)

max H0

(0)

H0

(0)

Heq

(0)

H0 (0)
x0 (t× 1 )

x0

x0 (t× 4 )

(0)

H1

(0)

Heq

(0)

H1 (t× 1 )

x1 (t× 1 )

x1

x1 (t× 4 )

Figure 10: Trajectoryn parameterisation
scenario

o example for a symmetric quadruple of spikes and transition 
 of five intervals
3
1
19
, 10
, 1 with cross-over moments t× = {8, 10, 20, 30}. Initial position is x0 , x1 = − 20
, − 34 .
characterised by γ = 1, 41 , 20




(0)
(0)
Equilibrium position is x0 , x1 =
− 43 , − 14 . The two rightmost spikes satisfy x2 (τ ) ≡ −x1 (τ ), H2 (τ ) ≡ H1 (τ ),
(0)

(0)

x3 (τ ) ≡ −x0 (τ ), H3 (τ ) ≡ H0 (τ ). Equilibrium height, equal for all spikes, is marked by the dashed line. Dotted lines show
1
the trajectory for γ ≡ 1. Diamonds mark the cross-over moments t× j . Other parameters used: (p, q, m, s) = (3, 2, 2, 0), ǫ = 10
.
Qualitatively similar results obtained with (p, q, m, s) = (2, 1, 2, 0), (2, 1, 3, 0), (4, 2, 2, 0).

28

(0)

max H0

(0)

H0

(0)

Heq

(0)

H0 (t× 1 )

x0 (t× 1 )

x0

x0 (t× 4 )

(0)

H1

(0)

Heq

(0)

H1 (t× 1 )

x1 (t× 1 )

x1

x1 (t× 4 )

29

(0)

H2

(0)

Heq

(0)

min H2

x2 (t× 4 )

x2

x2 (t× 1 )

(0)

H3

(0)

H3 (t× 1 )

(0)

min H3

x3 (t× 4 )

x3

x3 (t× 1 )

Figure 11:
example for a quadruple of spikes and transitionscenario of five
characterised

o
 intervals

n Trajectory parameterisation
3
1
19
9
.
, 10
, 1 with cross-over moments t× = {1, 2, 6, 8}. Initial position is x0 , x1 , x2 , x3 = − 20
, − 12 , 53 , 10
by γ = 1, 41 , 20

 

3
1 1 3
Equilibrium position is x0 , x1 , x2 , x3 = − 4 , − 4 , 4 , 4 . Equilibrium height, equal for all spikes, is marked by the dashed
line. Dotted lines show the trajectory for γ ≡ 1. Diamonds mark the cross-over moments t× j . Other parameters used:
1
(p, q, m, s) = (3, 2, 2, 0), ǫ = 10
. Qualitatively similar results obtained with (p, q, m, s) = (2, 1, 2, 0), (4, 2, 2, 0). With the set
(p, q, m, s) = (2, 1, 3, 0) the pattern failed to equilibrate.

H0

(0)

30

(0)

Heq

(0)

H0 (0)

19
− 20

−0.82
x0

(0)

H1

(0)

Heq

(0)

H1 (0)
− 21

−0.32
x1

31
(0)

H2

(0)

H2 (0)

(0)

min H2

0.32
x2

3
5

(0)

H3

(0)

H3 (0)

(0)

Heq

0.82
x3

9
10

Figure 12: Trajectory parameterisation example for a quadruple of spikes and inner layer transition scenario (figure 1) with




3
3
9
γ0 = , γ1 =
and ε = 5 in (6). Initial position is x0 , x1 , x2 , x3 =
− 19
. Equilibrium position is
, − 12 , 35 , 10
20
4


 20
x0 , x1 , x2 , x3 = − 34 , − 41 , 14 , 43 . Diamonds mark the cross-over moments t× j . Other parameters used: (p, q, m, s) =
1
(3, 2, 2, 0), ǫ = 10
. Qualitatively similar results obtained with (p, q, m, s) = (2, 1, 2, 0), (4, 2, 2, 0).

32

x′i (τ )

1

0

-1

τtrm

0
τ

(0)

Hi

(0)

max Hi

(0)

min Hi

-1

0
xi

Figure 13: Example of trajectory termination due to a bifurcation for a quadruple of spikes with γ ≡ 1. Top panel: velocity
x′i (τ ), i = 0 (leftmost spike, solid curve), i = 1 (left central spike, dashed curve), i = 2 (right central spike, dash-dotted
(0)

(0)

curve) and i = 3 (rightmost spike, dotted curve). Bottom panel: Hi (xi ) for the left spikes, i = {0, 1}, and Hi (−xi ) for
the
 right ones, i ={2, 3}. Diamonds
 mark trajectory termini.
 Curve linestyles
 are identical to top panel. Initial position is
9
− 19
. Equilibrium position x0 , x1 , x2 , x3 =
, − 12 , 53 , 10
20

− 43 , − 41 , 14 , 43 is never reached, the terminus
1
′
point τtrm characterised by infinite slope of xi (τ ). Other parameters used: (p, q, m, s) = (3, 3, 3, 0), ǫ = 10
.
x0 , x1 , x2 , x3 =

33
x′0 (τtrm )

x′1 (τ )

x′0 (τ )

x′1 (τtrm )

x′0 (0)

x′1 (0)

τtrm

τ× 1

τtrm

τ× 1

τ

τ

x′3 (τtrm )
x′3 (τ )

x′2 (τ )

x′2 (τtrm )

x′2 (0)

x′3 (0)
τtrm

τ× 1

τtrm

τ× 1

τ

τ

(0)

H0 (ttrm )
(0)

(0)

H1

H0

(0)

H1 (ttrm )

(0)

H0 (t× 1 )

(0)

H1 (t× 1 )
x0 (t× 1 )

x1

x0 (ttrm )

(0)

H2 (ttrm )

x1 (t× 1 )

x1 (ttrm )

x2

(0)

(0)

H3

H2

(0)

H3 (ttrm )

(0)

(0)

H3 (t× 1 )

H2 (t× 1 )
x2 (t× 1 )

x3

x2 (ttrm )

x4

x3 (t× 1 )

Figure 14: Example of trajectory termination due to a bifurcation for a quadruple of spikes with γ =

x3 (ttrm )

n

1, 14

o

and cross-over

moment t× = 1. Top: velocity x′i (τ ). Solid curves: i = 0 (leftmost spike), i = 1 (left central spike), i = 2 (right central spike)
(0)

and i = 3 (rightmost spike). Dotted curves show the same evolution for γ ≡ 1. Bottom: respective Hi (xi ). Diamonds mark

 

19
9
trajectory termini. Crosses mark the cross-over moment. Initial position is x0 , x1 , x2 , x3 = − 20
. Equilibrium
, − 12 , 53 , 10




3
1 1 3
position x0 , x1 , x2 , x3 = − 4 , − 4 , 4 , 4 is never reached, the terminus point τtrm characterised by infinite slope of x′i (τ )
1
.
for the the left pair of spikes, x′i (τtrm ) = 0 for the right pair. Other parameters used: (p, q, m, s) = (2, 1, 3, 0), ǫ = 10