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