Numerical solutions of a one-dimensional sub-diffusive Gierer-Meinhardt model with
controlled precision
Nehemie Nguimbous 1

Juan Pablo Ortiz Gil2

1Thompson Rivers University

Yana Nec 1

2Universidad de Antioquia

Subdiffusive Gierer-Meinhardt Model

Second Order Problem

This study focuses on sub-diffusion with Gierer-Meinhardt reaction kinetics. The GiererMeinhardt model is a mathematical model that describes the formation of spatial patterns in
biological systems. The sub-diffusive Gierer-Meinhardt model is given by

Once collecting terms of second order in equations (4a) and (4b) respectively, and upon postprocessing these equations, we obtain the following differential-algebraic system:

p
a
γ
∂t a = 2γ axx − a + q ,
h

m
a
γ
τ ∂t h = Dhxx − h + −γ s ,
h

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

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


−1 < x < 1,

t > 0,

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

(1)
(2)

dxi γ
dxi
−
sign
dσ
dσ





1. Establish a correspondence between the number of subdivisions n and the integral I
n(t,p,γ) → I(t,p,γ)(n).



I is computed using the composite Simpson method:

n
X


qbm 

B
m−s
B
m−s
=
H̄i
Gx(x; xj ) + H̄
hGxii f (p; γ),


(p + 1)H̄i  j=1

(7a)

H̄i(σ) = h(0)(x, t) = bm

n
X

h X
f (x)dx ≈
3
a

Z ∞
f (p; γ) =

H̄iB m−sG(x; xi),

Z ∞
bm =

umdy,



−1 < x < 1,

Gx(±1; xi) = 0,

up+1dy

(7b)

−∞

 . Z ∞

−∞

Solving the System: Matched Asymptotic Expansion
In order to fully capture the behavior of the activator and inhibitor, we scaled both the time
and position variables, respectively from
x − xi(σ)
γ+1
t to σ =  t, and x to yi =
,
γ

where the xis represent the positions where the activator and the inhibitor significantly interact. Matched asymptotic expansion involves the following steps:

du γ
Dy udy .
−∞ dy

(7c)

with D being a constant and


1 − −
+; x ) ,
hGxii = Gx (xi ; xi) + G+
(x
i
x i
2
γ
Dy u(y) = sign

Approximate the concentrations a and h with power series:



dxi
dσ



f (x2i) + 4f (x2i+1) + f (x2i+2) ,

b−a
h=
.
n

2. Our goal being to compute I with a controlled precision, it makes sense to establish a
correspondence between the number of subdivisions n and the residual R as well:
n(t,p,γ) → R(t,p,γ)(n).

Here, G is the solution of
DGxx − G = −δ(x − xi),

i

i=0

i=1

m .
p > 1, q > 0, m > 0, s ≥ 0, p−1
<
q
s+1

n/2−1 h

Z b

j6=i

where a(x, t) and h(x, t) are, respectively, the concentrations of the activator and inhibitor at
position x and time t. Here, 2γ and D denote the constant diffusivities, τ is the reaction time
constant, γ a real number such that 0 ≤ γ ≤ 1 and the exponents (p, q, m, s) satisfy

(0)
(1)
(2)
γ
−α
γ
2γ
A(yi, σ) = a(xi +  yi,  σ) = Ai (yi, σ) +  Ai (yi, σ) +  Ai (yi, σ) + · · ·

Procedure: Finding the Number of Subdivisions Required to Compute I
With a Precision of 10−10

This calculation is for (t, p, γ) values within the ranges:
{0.1, 1, 5} × {1.5, 2, 2.5, . . . , 4.5} × {0.1, 0.2, 0.3, . . . , 0.9}.
3. However, since the residual function we obtained is not a well-known function, we fit it
using a variant of the hyperbolic-arctangent function, denoted f = c1 × atanh(n) + c2, with
c1 and c2 being parameters used to minimize the distance between the residual and the
fitting curve. Consequently,

(8a)


    γ+1
Z ∞
1
dxi
1
u(y) − u y + sign
ξ1
dξ1,
Γ(−γ) 0
dσ
ξ1

n(t,p,γ) → f(t,p,γ)(n) ∼ R(t,p,γ).

(8b)

(3a)

Solutions
(0)

(1)

(2)

H(yi, σ) = h(xi + γ yi, −ασ) = Hi (yi, σ) + γ Hi (yi, σ) + 2γ Hi (yi, σ) + · · ·
Substitute these approximations into the original system and truncate the resulting
equations at a specific order (second order in our case):
γ

(0)

(1)

αγ ∂σ Ai + γ Ai

!

(0)

(1)

∼ ∂y2i Ai + γ Ai

!
−

(0)

(1)

!

Ai + γ Ai

(0)p
Ai
+
(0)q
Hi

(0)p
(1)
A
A
(1)
− qγ  i q 2 Hi + pγ i q .
(0)
(0)
Hi
Hi






(0)
(1)
(0)
(1)
(0)
(1)
γ
τ αγ ∂σ Hi γ Hi
∼ −2γ D∂y2i Hi + γ Hi
Hi + γ Hi



m
(0) 
(1)
(1) 
A
A
H
+ −γ i s 1 + γ m i − s i s  .
(0) 
(0)
(0)

H
A
H
i

i

(3b)

Upon solving the first equation of (7b) at x = xi, and for i = 1 (first spike), we have

 −1
B m−1−s
H̄1 = bmG(x1; x1)
.
(0)

(4a)

where n(t,p,γ) and R(t,p,γ) represent the number of subdivisions and the corresponding
residual values for specific t, p, and γ values.
(4b)

Verification of the Accuracy

i

The table below displays the approximated number of subdivisions obtained using the inverse of
the fitting curve for t = 0.1 and some discrete values of p and γ.
Figure 1. Activator concentration for p = 2.

For α = γ + 1, and upon collecting terms of first order in equations (4a) and (4b), we respectively
obtain:
(0)p
Ai
(0)
(0)
2
2 H (0) = 0,
∂yi Ai (yi) − Ai +
=
0,
∂
−∞ < yi < ∞.
(5)
yi i
(0)q
Hi
Upon solving this system, one obtains
q
with g being a function of the variable σ, B =
and
p−1


u(y) =

 1
p+1
(p − 1)y p − 1
2
sech
.
2
2

Figure 2. Inhibitor concentration for p = 2.

Computing Dtγ u

First Order Problem

(0)
Ai (σ) = H̄iB(σ)u(yi),

4. To find the number of subdivisions based on a specific residual (the inverse path), we
determine the inverse f −1 of f such that
−1 (R) → n
f(t,p,γ)
(t,p,γ),

Collect terms of the same order: the first and second order terms are collected separately
to form the first and second order problems.

(0)
Hi (yi) = H̄i(σ) = g(σ),

(0)

Substituting this result into the second equation in (6) and computing both A1 and H1 = H̄1
lead to the following solutions:

Figure 3. Specific case for t = 1, p = 4, and γ = 0.8. The x-axis represents the number of subdivisions and the
y-axis the logarithm of the corresponding residuals. The solid curve represents the numerical results while
the dashed one represents the fitting curve.

(6)

dxi
γ
In order to compute
and track the evolution of xi in (7a), we have to compute Dt u. We first
dσ
start by eliminating the singularity at ξ1 = 0 by regularizing its expression. This leads to
−γ
t∞

1−γ
t∞

1
γ
0
Dt u(t) =
u(t) − u(t − t∞) +
u (t − t∞) −
Γ(1 − γ)
Γ(2 − γ)
Γ(2 − γ)





Z t−t∞

u00(y)(t − y)1−γ dy.

t

γ
The first two terms of Dt u are known values, except for the third one, named I, which is the

integral of a complex function and is impossible to compute analytically. Moreover, I depends
on three parameters: p, γ, and t. The number of subdivisions required to compute I varies with
these parameters as well. We aim to predict the number of subdivisions n required to compute
I with a precision of 10−10.

g/p
1.5
2.0
2.5
3
3.5
4
4.5
0.1 35065
51414
61242
73570
80108
86651
93201
0.3 214565
342453
430919
519457
563799
608143
686814
0.5 1674913 2341552
3365120
3698389
4365537
5033725
5366720
0.7 4662043 28071526 37953058
45866233
54342976
60260601
68207443
0.9 76155108 727802772 1015373018 1025834015 1025841827 1025918580 1025859852

Table 1. Approximated number of subdivisions for t = 0.1

As a test for the effectiveness of our method, we computed the absolute value of the difference
γ
between Dt u

N =2n

γ
and Dt u

for some values of t. This leads to the following table:
N =n

t
3.6
3.7
3.8
3.9
4.0
error -1.784137e-010 -1.721758e-010 -1.664409e-010 -1.605218e-010 -1.547047e-010
t
4.1
4.2
4.3
4.4
4.5
error -1.494032e-010 -1.443059e-010 -1.395319e-010 -1.351115e-010 -1.307046e-010
Table 2. Error file for p = 2 and γ = 0.7.

CAIMS 2024, Kingston ON

nguimbousn21@mytru.ca