A Sufficient Condition for Asymptotic Stability in a Modified Fowler Model and Seashell Pattern Reconstruction Using PINNs

A Sufficient Condition for Asymptotic Stability in a Modified Fowler Model and Seashell Pattern Reconstruction Using PINNs Jingkai Wang∗ Abstract We investigate seashell pattern formation through an activator-inhibitor system, focusing particularly on a modified Fowler model. The modification in the activator dynamics preserves the biological mechanisms while enhancing the tractability for stability analysis. The finite element scheme is developed to solve the system numerically on irregular seashell surfaces. We propose two key theorems: the first demonstrates a sufficient condition on the parameters for the existence of a stable fixed point using Turing instability analysis; the second derives an amplification matrix that characterizes the stability region of the finite element scheme. The parameter conditions for the existence of a Hopf bifurcation in the modified Fowler model without diffusion are derived and verified. Given parameters satisfying the sufficient condition for asymptotic stability, numerical simulations confirm the convergence behavior of the dynamical system over time. We explore the inverse parameter identification using physics-informed neural networks (PINNs), capturing the overall seashell patterns effectively. Keywords: seashell patterns, Fowler model, finite element method, variational formulation, Turing instability analysis, Hopf bifurcation, stability region, physics-informed neural networks 1 Introduction The pigmentation of seashells manifests complex and heterogeneous spatial patterns, deriving from the biochemical effects along the growth edge. To understand the development of seashell patterns, reaction-diffusion models such as activator–inhibitor systems have been proposed [7]. The activator-inhibitor system model introduces two interacting dynamical systems: the first captures the activator that facilitates the growth of spatial patterns, while the second system denotes the inhibitor that suppresses the activator development. Two activator-inhibitor system models employed in seashell pattern studies are the Gierer–Meinhardt model and a modified Fowler model. Comparing to the Fowler model [6], the modified version in our work reconstructs the nonlinear term of the activator dynamics to preserve biological meanings and improve the tractability for stability analysis. While the Gierer–Meinhardt model is effective in capturing the activator–inhibitor dynamics in seashell patterns, the modified Fowler model achieves more realistic pattern simulations by introducing the global regulatory effects. A wide range of studies has applied the finite element methods (FEM) and other numerical methods on the reaction-diffusion model [3, 10, 26]. The stability condition for the reduced Gierer–Meinhardt (GM) model has been introduced, where the fixed point can be derived analytically [2]. The physics-informed neural networks have been proposed as a framework for solving inverse problems, including the identification of parameters and initial conditions of partial differential equations (PDE) [21]. ∗Department of Applied Physics and Applied Mathematics, Columbia University ([email protected]) 1 In our work, we leverage finite element methods to derive the numerical results of the 2-dimensional activator-inhibitor systems. After the variational formulation, the non-linear system is discretized using a time-stepping scheme such as a semi-implicit scheme [17, 18]. Since the activator-inhibitor system captures the interaction behaviors, the numerical solution is updated at each time step until convergence. Two key theorems are proposed. Theorem 1 is derived from the Turing instability analysis, yielding a sufficient condition to analyze the convergence behavior of the activator-inhibitor system. Although the fixed point cannot be computed in closed form and may not exist in general, the sufficient condition for its existence and stability is derived analytically. Theorem 4 derived an amplification matrix based on Theorem 1 to capture the stability region of the finite element scheme. We have rigorously analyzed the Hopf bifurcation conditions by deriving the parameter threshold and associated inequalities in Theorem 2 and Theorem 3. We implement the numerical methods on the irregular seashell shape using the parameters that satisfy the sufficient condition. The physics-informed neural networks (PINNs) are introduced to identify the parameters from the real seashell image. Figure 1: A image of seashell patterns. 2 Activator Inhibitor Systems for Seashell Pattern Modeling As a subclass of reaction-diffusion models, the activator-inhibitor system is leveraged to measure the evolution of spatial patterns over time [11]. A real image of seashell patterns is displayed in Figure 1. Here, u(x, y, t) denotes the density of the activator that promotes pattern formation, while v(x, y, t) denotes the inhibitor concentration responsible for suppressing the activator population. At each time step, the distribution of uis employed to represent the pigment intensity reflecting the surface patterns of seashells. The Gierer-Meinhardt model and the modified Fowler model are introduced and analyzed. 2.1 Model Utilization The Giere-Meinhardt (GM) Model [7] measures the biomedical interaction between activator and inhibitor using the following formula: ∂tu=Du∇2u+ρu2 v(1 + κu2)−µu,(1) ∂tv=Dv∇2v+ρ(u2−kuv),(2) 2 where Du, Dv, ρ, κ, µ > 0. In the formula, parameter ρcaptures the speed of reaction kinetics between activator and inhibitor; κrepresents the saturation constant that modulates the growth of the activator; µis the activator decay rate and kuis the inhibitor decay rate. In the diffusion term, parameters Duand Dvrepresent spatial diffusion rates of the activator and inhibitor respectively. Homogeneous solutions appear when the activator diffuses rapidly (large Du); while the imbalance pattern formulation due to Turing instability becomes prominent when the inhibitor exhibits a higher diffusion rate (large Dv). The emergence of biological patterns is determined by the parameters that reflect the reaction, diffusion and interaction behaviors, along with the initial state. The Fowler model [6] is motivated by the Giere-Meinhardt model by retaining activator-inhibitor dynamics, while the activator-inhibitor system model leverages more parameters to prevent unbounded growth and enhance robustness. In this work, we introduce a modified Fowler model where the reaction term in the activator dynamics is reformulated. The modification preserves the biological mechanisms of the activator-inhibitor system while suppressing the strength of the non-linear feedback. Among new parameters, ρ0controls the activator growth by dampening the nonlinear feedback, σprevents the inhibitor concentration from collapsing to zero level; h0is leveraged to prevent singular behavior. The modified Fowler model is written as ∂tu=Du∇2u+ρ (v+h0)(1 + κu2+ρ0)u2−µ2u, (3) ∂tv=Dv∇2v+σ+ρu2 1 + κu2−ηv, (4) where Du, Dv, µ2, η, σ, ρ, ρ0, κ, h0>0. Both the Giere-Meinhardt model and the modified Fowler model are used to capture the variation of pigmentation patterns of seashells at each time step by measuring the biological interaction between activator and inhibitor. In certain boundaries, the activator and inhibitor demonstrate time-dependent dynamics, formulated by the equations, providing the spatial patterns. The finite element method (FEM) with triangular discretizations is leveraged to reduce the IBVP problem into a system of ordinary differential equations (ODEs), which can be evolved in time using implicit schemes. Specifically, when parameters ρ, σ, h0are set to be zero, the modified Fowler model is reduced to the Giere-Meinhardt model. 2.2 A Sufficient Condition for Stable Fixed Point As shown in [2], a stability condition was derived for a reduced Gierer-Meinhardt model. However, the modified Fowler model is more intricate since the fixed point cannot be computed analytically. To analyze the stability of fixed points in the activator-inhibitor system of the modified Fowler model, we consider the following system without the diffusion term. ∂tu=ρ u2 (v+h0) (1 + κu2+ρ0)−µ2u, (5) ∂tv=σ+ρu2 1 + κu2−ηv. (6) The fixed point (u∗, v∗) for the above non-linear system satisfies the following equations that ρu∗=µ2(v∗+h0)1 + κu∗2+ρ0(7) v∗=1 ησ+ρu∗2 1 + κu∗2(8) 3 By replacing the formula of v∗, it can be computed that u∗should satisfy the equation that f(u∗) = µ2 u∗1 ησ+ρ·u∗2 1 + κu∗2+h0·1 + κu∗2+ρ0−ρ= 0.(9) And the Jacobian matrix at the fixed point can be derived as J(u∗, v∗) = a11 a12 a21 a22=    µ22(1 + ρ0) 1 + κ(u∗)2+ρ0 −1−µ2u∗ v∗+h0 2ρu∗ (1 + κ(u∗)2)2−η     (10) To study the Turing instability, we analyze the behaviors of spatial perturbations evolving around the steady state at (u∗, v∗). Suppose Dis the diffusion matrix and d=Du/Dv, and then the Jacobian matrix is modified by incorporating the effect of diffusion with Jk=J−k2D: Jk=    µ22(1 + ρ0) 1 + κ(u∗)2+ρ0 −1−dk2−µ2u∗ v∗+h0 2ρu∗ (1 + κ(u∗)2)2−η−k2     , D =d0 0 1.(11) Since the stability of the fixed points requires the real part of all eigenvalues to be negative, this condition is equivalent to ensuring that Tr(Jk)<0 and Det(Jk)>0 for non-negative integer k. Tr(Jk) = µ22(1 + ρ0) 1 + κ(u∗)2+ρ0 −1−η−(d+ 1)k2<0,(12) Det(Jk) = µ22(1 + ρ0) 1 + κ(u∗)2+ρ0 −1−dk2·(−η−k2) + µ2u∗ v∗+h0 ·2ρu∗ (1 + κ(u∗)2)2>0.(13) One sufficient condition to ensure the inequalities is requiring u∗>q1+ρ0 κ. From Equation (9), we observe that f(u∗)→ ∞ as u∗→ ∞. Therefore, if f(q1+ρ0 κ)<0, by the Intermediate Value Theorem, there exists at least one fixed point satisfying u∗>q1+ρ0 κ. Then Theorem 1 can be derived to demonstrate the sufficient condition. Theorem 1. There exists a fixed point (u∗, v∗), and the fixed point is stable if the following condition holds, where 0< µ2< ρ·r1 + ρ0 κ (2 + 2ρ0)1 ησ+ρ(1 + ρ0) κ(2 + ρ0)+h0.(14) Moreover, if the sufficient condition holds, there exists a fixed point that is asymptotically stable. In the small neighborhood of the fixed point (u∗, v∗), the activator concentration and inhibitor concentration converge to it over time. 4 2.3 Hopf Bifurcation Analysis The occurrence of Hopf bifurcation in reaction-diffusion systems is widely analyzed [15, 23]. The Hopf bifurcation condition for the modified Fowler model is analyzed based on the reaction terms of the system. Then the Hopf bifurcation occurs when a pair of complex conjugate eigenvalues cross the imaginary axis, i.e. Tr(J0) = 0 and Det(J0)>0. From Equation (12), the appearance of Hopf bifurcation requires µ2> η, and the steady state is u∗= ˆu∗=s(µ2−η)(ρ0+ 1) κ(µ2+η).(15) Since the steady state of the system satisfies Equation (9), we can derive a parameter condition where the system undergoes a Hopf bifurcation. Theorem 2. If the Hopf bifurcation occurs in the modified Fowler model at the fixed point (ˆu∗,ˆv∗), then the parameters satisfy the formula that h0=ˆ h0=ρ 2µ2 2(ρ0+ 1)r(ρ0+ 1)(µ2−η)(η+µ2) κ−1 ησ+ρ(ρ0+ 1)(µ2−η) κ(ρ0(µ2−η)+2µ2).(16) Derived from Equation (9), Equation (16) enforces Tr(J0)=0. However, the occurrence of a Hopf bifurcation additionally requires that Det(J0)>0,∆ = Tr(J0)2−4Det(J0)<0, and ∂ ∂h0Tr(J0)= 0. The derivative of Tr(J0) in terms of h0cannot be computed directly, since Tr(J0) depends implicitly on u∗, which satisfies Equation (9). To determine whether ∂ ∂h0Tr(J0)= 0, we analyze how changes in h0influence the steady state u∗. Corollary 1. In the non-linear dynamical system corresponding to the modified Fowler model without diffusion terms, the condition ∂ ∂h0Tr(J0)= 0 is equivalent to the condition ∂f ∂x (u∗, h0)= 0, ensuring the transversality condition for a Hopf bifurcation. Proof. From Equation (9), we extend the function fto f1(x, y), s.t. f1(u∗, h0) = f(u∗) = 0. The function f1is defined as f1(x, y) = µ2 x1 ησ+ρ·x2 1 + κx2+y·1 + κx2+ρ0−ρ= 0.(17) Then f1is continuously differentiable in the neighborhood of (u∗, h0) and ∂f1 ∂x (u∗, h0)= 0. According to the implicit function theorem (IFT) [12], there exists a continuously differentiable function u∗(h0) s.t. f1(u∗(h0), h0) = 0. Then it can be derived that ∂ ∂h0 Tr(J0) = −∂ ∂u∗Tr(J0)∂f1(u∗, h0)/∂y ∂f1(u∗, h0)/∂x.(18) Since ∂ ∂u∗Tr(J0)= 0, the condition ∂ ∂h0Tr(J0)= 0 holds if and only if ∂f ∂x (u∗, h0)= 0. From Theorem 2 and Corollary 1, the remaining conditions for the existence of a Hopf bifurcation in the modified Fowler model can be derived. 5 Theorem 3. Suppose the steady state of the reaction system satisfies Equation (15), parameters fulfill the condition in Equation (16) and µ2> η. A Hopf bifurcation arises if the following criteria are met: (i) The Jacobian determinant Det(J0)>0at the steady state is positive, i.e., r−(η2−µ2 2)(1 + ρ0) κ(ηρ0−µ2(2 + ρ0))2−4µ5 2(1 + ρ0)2+η4ρ2 0qκ(−η2+µ2 2)(1 + ρ0) −2η3µ2ρ0qκ(−η2+µ2 2)(1 + ρ0)(2 + ρ0)2+η2µ2 24µ2(1 + ρ0)2+qκ(−η2+µ2 2)(1 + ρ0)(2 + ρ0)2<0. (19) (ii) The eigenvalues are complex when the discriminant is negative, requiring ∆ = Tr(J0)2− 4Det(J0)<0, i.e., µ22(ρ0+ 1)(η+µ2) 2µ2(ρ0+ 1) −1−η2 −4  2µ3 2ρ(ρ0+ 1)2(µ2−η)(η+µ2) κρ (η+µ2+ (ρ0+ 1)(µ2−η))2q(ρ0+1)(µ2−η)(η+µ2) κ −ηµ22(ρ0+ 1)(η+µ2) 2µ2(ρ0+ 1) −1<0. (20) (iii) The derivative of the trace ∂ ∂h0Tr(J0)= 0 is non-zero, which is equivalent to ∂f ∂x (u∗, h0)= 0, i.e., −4µ5 2(1 + ρ0)2+η4ρ2 0qκ(−η2+µ2 2)(1 + ρ0)−2η3µ2ρ0qκ(−η2+µ2 2)(1 + ρ0)(2 + ρ0) +η2µ2 24µ2(1 + ρ0)2+qκ(−η2+µ2 2)(1 + ρ0)(2 + ρ0)2= 0.(21) To rigorously observe a Hopf bifurcation in the system, the criteria in Theorems 2 and 3 must be strictly satisfied. We consider the range defined by hmin =ˆ h0−ϵ1and hmax =ˆ h0+ϵ2where ϵ1, ϵ2>0. Within the parameter range, a Hopf bifurcation might occur in the system. However, for certain values h0, Equation (9) might admit no real solutions, indicating no fixed point in the system. In such cases, the conditions for the Hopf bifurcation are not well defined. The existence of ϵ1and ϵ2is verified in Corollary 2. Corollary 2. In this system, there exist ϵ1, ϵ2>0such that a Hopf bifurcation is observed when h0∈(hmin, hmax)where hmin =ˆ h0−ϵ1and hmax =ˆ h0+ϵ2. Note that ˆ h0is defined in Equation (16). Proof. From Equation (17), it is equivalent to show that ∀h0∈(hmin, hmax), ∃u∗>0 such that f1(u∗, h0) = 0. According to Corollary 1 and the implicit function theorem (IFT), there exists a unique function such that g(ˆ h0) = ˆu∗and f1(g(h0), h0) = 0 in the neighborhood of (ˆu0,ˆ h0). From Figures 2 and 3, we vary the parameter h0within a small neighborhood of the critical value ˆ h0, where ϵ1=ϵ2= 0.0001 are employed to see different behaviors. When h0<ˆ h0, the fixed point is a stable spiral; when h0>ˆ h0, the fixed point is an unstable spiral and a stable limit cycle emerges. The behaviors of trajectories indicate the occurrence of a supercritical Hopf bifurcation. 6 (a) Stable fixed point when h0=ˆ h0− 0.0001. (b) Distance to fixed point in different trajectories. Figure 2: Hopf bifurcation observation when h0=ˆ h0−0.0001. Parameters: µ2= 0.02, η= 0.0155, σ= 0.001, ρ= 0.2, ρ0= 0.1, κ= 0.2847, and h0= 0.1, with time step ∆t= 0.01 (a) Unstable fixed point and stable limit cycle when h0=ˆ h0+ 0.0001. (b) Distance to fixed point in different trajectories. Figure 3: Hopf bifurcation observation when h0=ˆ h0+ 0.0001. Parameters: µ2= 0.02, η= 0.0155, σ= 0.001, ρ= 0.2, ρ0= 0.1, κ= 0.2847, and h0= 0.1, with time step ∆t= 0.01. 3 Finite Element Formulation To determine the numerical solutions of the activator-inhibitor system, the finite element method (FEM) is employed. The variational formulation is used to derive the weak form by multiplying each equation with a test function, as shown in Equations (23) and (24). Then the spatial domain is discretized using triangular meshes, as displayed in Figure 4. The system of partial differential equations (PDEs) is posed as an initial boundary value problem (IBVP), where the initial condition is prescribed at time t= 0. To ensure biological significance, we assume that no biochemical instances flow across the domain boundary. Since the system is constrained to the seashell surface, we adopt the Neumann boundary condition in terms of the outward normal vector n, as commonly applied in pattern formulation and reaction-diffusion models [9, 20, 19]: ∂u ∂n =∂v ∂n = 0,on ∂Ω.(22) Among the two systems, the modified Fowler model extends the Giere-Meinhardt (GM) system by incorporating additional parameters. Therefore, we derive the numerical methods for the modified 7 Fowler model, which simplifies to the GM model when certain parameters vanish. 3.1 Variational Formulation and Time Discretization To consider the variational formulation of the modified Fowler model, we multiply each PDE with a test function and integrate both sides over the domain. After using Green’s 1st identity, the zero-flux conditions are leveraged to derive the weak form, where the boundary terms vanish due to homogeneous Neumann boundary conditions [16]: ZΩ ∂tu ϕ dx +DuZΩ ∇u· ∇ϕ dx =ZΩρu2 (v+h0)(1 + κu2+ρ0)−µ2uϕ dx, (23) ZΩ ∂tv ψ dx +DvZΩ ∇v· ∇ψ dx =ZΩσ+ρu2 1 + κu2−ηvψ dx. (24) The test functions ϕand ψbelong to the Sobolev space H1(Ω) [5, 16], suggesting that the first derivatives of ϕand ψare square integrable. The use of test functions enables the extraction of residual dynamics of the weak form. Since the test functions are defined on the same Sobolev space, we leverage the same basis {φi(x, y)}m i=1 for the finite element subspace Vh⊂H1(Ω) [13]. Functions u(x, y, t) and v(x, y, t) can be approximated in the finite element space [16, 24]: uh(x, y, t) = m X i=1 Ui(t)φi(x, y), vh(x, y, t) = m X i=1 Vi(t)φi(x, y).(25) Several authors apply Galerkin finite element formulations to reaction-diffusion models [1, 14]. By enforcing the variational form on the basis function φj∈Vh, the corresponding finite element equations in the vector form can be derived as MdU dt +DuKU =F(U, V ),(26) MdV dt +DvKV =G(U, V ),(27) where we define the mass matrix M, the stiffness matrix Fand non-linear vectors F(U, V ) and G(U, V ). The mass matrix Mcaptures the inner product of finite element basis; the stiffness matrix Fcorresponds to the diffusion term using the gradients of basis functions; the non-linear vectors F(U, V ) and G(U, V ) are derived from the reaction terms of the modified Fowler model, approximated by uhand vh. Mij =ZΩ φiφjdx, Kij =ZΩ ∇φi· ∇φjdx, (28) Fj(U, V ) = ZΩρu2 h (vh+h0)(1 + κu2 h+ρ0)−µ2uhφjdx, (29) Gj(U, V ) = ZΩσ+ρu2 h 1 + κu2 h −ηvhφjdx. (30) The above ODE systems involve the first-order temporal derivatives and non-linear interactions between variables. To further solve this equation numerically, authors use a time discretization scheme [27]. Despite the computational simplification, the forward Euler method requires extremely 8 small time steps to ensure stability due to the stiffness of the system. The backward Euler method is adopted to support stable simulations over larger time intervals [8]. MUn+1 −Un ∆t+DuKUn+1 =F(Un+1, V n+1),(31) MVn+1 −Vn ∆t+DvKV n+1 =G(Un+1, V n+1).(32) Suppose that uh(x, y, tn) can be approximated by Pm i=1 Un iφi(x, y) at time t=tn, then Unrepresents the coefficients of the basis functions, similar to Vn. Un= [Un 1, Un 2,· · · , Un m]T, V n= [Vn 1, V n 2,· · · , V n m]T.(33) We assume Au=M+ ∆tDuK,Av=M+ ∆tDvK,bu=MUnand bv=MV n. To facilitate efficient computation, the authors leverage the semi-implicit scheme, where the nonlinear reaction terms are evaluated at the current time step [18]. AuUn+1 =bu+ ∆t F(Un, V n),(34) AvVn+1 =bv+ ∆t G(Un, V n).(35) Overall, the variational formulation converts the activator-inhibitor system into a weak formulation with test functions in H1(Ω). The solutions are expressed by the finite element basis functions that define the subspace of Sobolev space. Spatial discretization using the basis leads to a system of coupled vector ODE equations. We adopt the semi-implicit scheme to treat the diffusion term implicitly and the reaction term explicitly [17, 18]. The scheme maintains temporal stability and addresses the stiffness of the system. Alternatively, the method is interpreted as a backward Euler scheme combined with Picard linearization, where the non-linear algebraic system is approximated using Picard iterations by replacing F(Un+1, V n+1) and G(Un+1, V n+1) with F(Un, V n) and G(Un, V n), resolving the interaction system numerically. 3.2 Triangular Discretization Figure 4: Triangular mesh discretization for the seashell boundary with 500, 2000, and 10,000 mesh points, providing spatial resolution for FEM simulations. Since we are simulating the pigmentation patterns in seashells, the 2-dimensional shape of the 9 [14] Siqing Li, Leevan Ling, Steven J Ruuth, and Xuemeng Wang. Realistic pattern formations on surfaces by adding arbitrary roughness. SIAM Journal on Applied Mathematics, 84(3):1163– 1185, 2024. [15] Yunfeng Liu and Yuanxian Hui. Hopf bifurcation in a delayed reaction–diffusion–advection equation with ideal free dispersal. Boundary Value Problems, 2021:1–20, 2021. [16] Anders Logg, Kent-Andre Mardal, and Garth Wells. Automated solution of differential equations by the finite element method: The FEniCS book, volume 84. Springer Science & Business Media, 2012. [17] Anotida Madzvamuse. Time-stepping schemes for moving grid finite elements applied to reaction–diffusion systems on fixed and growing domains. Journal of computational physics, 214(1):239–263, 2006. [18] Anotida Madzvamuse. A modified backward euler scheme for advection-reaction-diffusion systems. Mathematical Modeling of Biological Systems, Volume I: Cellular Biophysics, Regulatory Networks, Development, Biomedicine, and Data Analysis, pages 183–189, 2007. [19] Philip K Maini and Thomas E Woolley. The turing model for biological pattern formation. The dynamics of biological systems, pages 189–204, 2019. [20] Takashi Miura and Philip K Maini. Periodic pattern formation in reaction—diffusion systems: An introduction for numerical simulation. Anatomical Science International, 79:112–123, 2004. [21] Maziar Raissi, Paris Perdikaris, and George E Karniadakis. Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations. Journal of Computational physics, 378:686–707, 2019. [22] Jos´e A Rodrigues and Beatriz Vieira. A mathematical analysis of image mesh generation using delaunay triangulation and image processing techniques. In 4th ROME International Conference on Challenges in Engineering, Medical, Economics and Education: Research & Solutions (CEMEERS-24b), 2024. [23] Yongli Song, Heping Jiang, and Yuan Yuan. Turing-hopf bifurcation in the reaction-diffusion system with delay and application to a diffusive predator-prey model. J. Appl. Anal. Comput, 9(3):1132–1164, 2019. [24] Endre S¨uli. Lecture notes on finite element methods for partial differential equations. Mathematical Institute, University of Oxford, 2012. [25] Vidar Thom´ee. Galerkin finite element methods for parabolic problems, volume 25. Springer Science & Business Media, 2007. [26] G Wu, Eric Wai Ming Lee, and Gao Li. Numerical solutions of the reaction-diffusion equation: An integral equation method using the variational iteration method. International Journal of Numerical Methods for Heat & Fluid Flow, 25(2):265–271, 2015. [27] Congcong Xie and Xianliang Hu. Finite element simulations with adaptively moving mesh for the reaction diffusion system. Numerical Mathematics: Theory, Methods and Applications, 9(4):686–704, 2016. 16