Turing Pattern Formation Under Heterogeneous Distributions of Parameters for an Activator-Depleted Reaction Model

Abstract

The purpose of this article is to study Turing pattern formation in one- and two-dimensional domains under heterogeneous distributions of the parameters for an activator-depleted model. Unlike previous studies of this nature, the choice of the heterogeneous distributions of the parameters is closely linked and estimated by use of rigorous wave mode selection in order to excite different modes in different subsets of the domains. This allows us to relate the numerical solutions with theoretical linear stability analytical results. Our most revealing results show that the wave modes of adjacent subsets evolve locally and yet possess continuity across the interface. These local patches of the solutions result in a globally heterogeneous solution stable only in the presence of heterogeneous distributions of the parameters. Furthermore, our results show that initial conditions continue to play a crucial role in the selection of excitable wave modes and consequently the formation of the inhomogeneous patterns formed. In particular, initial conditions influence pattern orientation and polarity, and yet with a prepattern, the patterns conserved orientation and polarity. Numerical solutions are obtained by the use of the finite element method and the backward Euler scheme to deal with the spatial and the time discretisations, respectively.

Appendices

Appendix A: Numerical Solution

1.1 Appendix A.1: The Finite Element Weak Formulation

Let $w\in H^1(\Omega )$ be a test function. Integrating over $\Omega$ the product of Eq. (1) and w, the weak formulation following integration by parts reads: find $u,v\in L^2\big ([0,T],H^1(\Omega )\big )$ such that:

$$\begin{aligned} \begin{aligned}&\begin{aligned} \int _\Omega w\frac{\partial u}{\partial t}d\Omega +\int _\Omega \nabla w\cdot \nabla u d\Omega -\int _\Omega w\gamma f(u,v)d\Omega -\int _\Gamma (\nabla u.{\mathbf {n}})w d\Gamma =0, \end{aligned} \\&\begin{aligned} \int _\Omega w\frac{\partial v}{\partial t}d\Omega +\int _\Omega \nabla w\cdot d\nabla vd\Omega -\int _\Omega w\gamma g(u,v)d\Omega -\int _\Gamma (\nabla v.{\mathbf {n}})w d\Gamma =0, \end{aligned} \end{aligned} \end{aligned}$$

(A.1)

where ${\mathbf {n}}$ is the normal vector to $\Gamma$. Bringing in the zero Neumann boundary conditions, the last term of Eq. (A.1) vanishes.

1.2 Appendix A.2: Spatial Discretisation

To build the finite element approximation, let us consider $\Omega ^h\subset \Omega$ a discretisation of $\Omega$ with vertices $\varvec{x}^h_j\in \Omega$, $j=1,\ldots ,N_n$ defined by: $\Omega ^h=\bigcup _{e=1}^{N_e}\Omega _e^h$. In one-dimensional domains, $\Omega _e^h$, $e=1,\ldots ,N_e$, are segments limited by two vertex $\varvec{x}^h_j$, while in two-dimensional domains are quadrilaterals defined by four vertex $\varvec{x}^h_j$. $\Omega _e^h$ are defined such that: $\bigcap _{e=1}^{N_e}\text {Int}\Omega _e^h=\emptyset$, where $\text {Int}\Omega _e^h$ is the interior of $\Omega _e^h$ and $\emptyset$ denotes the empty set.

Let us now define the finite element space as:

$$\begin{aligned} S^h(\Omega ^h)=\big \{\varphi \in&C^0(\Gamma ^h):\varphi \big |_{\Omega _e^h} \text { is linear affine for each }\Omega _e^h\in \Omega ^h\big \}, \end{aligned}$$

with basis $\{\chi _i\}$ $i=1,\ldots ,N_n$, such that $\chi _i(\varvec{x}_j^h)=\delta _{ij}$ ($\delta$ the Kronecker delta).

Consider the piece-wise linear approximations of $u(\cdot ,t)$ and $v(\cdot ,t)$, respectively, $u^h(\cdot ,t)\in S^h(\Omega ^h)$ and $v^h(\cdot ,t)\in S^h(\Omega ^h)$ defined by:

$$\begin{aligned} u^h(\cdot ,t)=\sum _{j=1}^{N_n}\chi _j(\cdot )U_j^h(t),\quad \text {and} \quad v^h(\cdot ,t)=\sum _{j=1}^{N_n}\chi _j(\cdot )V_j^h(t) \end{aligned}$$

(A.2)

and the general form of $\varphi$ as:

$$\begin{aligned} \varphi (\cdot ,t)=\sum _{i=1}^{N_n}\chi _i(\cdot )\eta _i(t), \end{aligned}$$

(A.3)

where $\eta _i(t)$ are arbitrary bounded values.

Then, the discrete form of the weak problem Eq. (A.1) is: find $u^h,v^h\in S^h(\Omega ^h)$ such that

$$\begin{aligned} {\left\{ \begin{array}{ll} \frac{\partial U_j}{\partial t}\int _{\Omega ^h}\chi _i\chi _j d{\Omega ^h} +U_j\int _{\Omega ^h} \nabla \chi _i\cdot \nabla \chi _j d{\Omega ^h} =\int _{\Omega ^h} \chi _i\gamma f(\chi _jU_j,\chi _jV_j)d{\Omega ^h}, \\ \\ \frac{\partial V_j}{\partial t}\int _{\Omega ^h} \chi _i\chi _jd{\Omega ^h} +V_j\int _{\Omega ^h} \nabla \chi _i\cdot d\nabla \chi _jd{\Omega ^h} =\int _{\Omega ^h} \chi _i\gamma g(\chi _jU_j,\chi _jV_j)d{\Omega ^h}, \end{array}\right. } \end{aligned}$$

(A.4)

notice that the sums for $u^h$, $v^h$ and $\varphi$, the spatial dependency of $\chi _i$ and $\chi _j$, and the time dependency of $U_j(t)$ and $V_j(t)$ are omitted for notation simplicity, and that since $\eta _i(t)$ are arbitrary they cancel out.

Let us rewrite Eq. (A.4) as:

$$\begin{aligned} \begin{aligned}&\sum _{e=1}^{N_e}\left( {\mathbb {M}}_e\cdot \frac{\partial {\mathbf {U}}_e}{\partial t}+{\mathbb {K}}_e\cdot {\mathbf {U}}_e-\gamma {\mathbf {F}}_e \right) ={\mathbf {0}},\\&\sum _{e=1}^{N_e}\left( {\mathbb {M}}_e\cdot \frac{\partial {\mathbf {V}}_e}{\partial t}+d{\mathbb {K}}_e\cdot {\mathbf {V}}_e-\gamma {\mathbf {G}}_e \right) ={\mathbf {0}}, \end{aligned} \end{aligned}$$

(A.5)

where ${\mathbf {U}}_e$ and ${\mathbf {V}}_e$ are local nodal values of $\Omega ^h_e$, ${\mathbb {M}}_e$ are local mass matrices:

$$\begin{aligned} {\mathbb {M}}_e=\int _{\Omega ^h_e} \chi _i\chi _jd{\Omega ^h_e}, \end{aligned}$$

(A.6)

${\mathbb {K}}_e$ are local stiffness matrices:

$$\begin{aligned} {\mathbb {K}}_e=\int _{\Omega ^h_e} \nabla \chi _i\cdot d\nabla \chi _jd{\Omega ^h_e}, \end{aligned}$$

(A.7)

and ${\mathbf {F}}_e$ and ${\mathbf {G}}_e$ are local reaction vectors:

$$\begin{aligned} {\mathbf {F}}_e=\int _{\Omega ^h_e} \chi _i f(\chi _jU_j,\chi _jV_j)d{\Omega ^h_e},\quad \text {and} \quad {\mathbf {G}}_e=\int _{\Omega ^h_e} \chi _i g(\chi _jU_j,\chi _jV_j)d{\Omega ^h_e}.\qquad \end{aligned}$$

(A.8)

1.3 Appendix A.3: The Temporal Discretisation Element-Wise

So far the system is completely discrete in space. Now, let us deal with the time-dependent terms. Consider $\tau >0$ a time step such that $t^k=k\tau$ with $k=1,\ldots ,N_T$ and $N_T$ the number of steps to discretise [0, T]. Applying the backward Euler scheme on Eq. (A.4) the fully discrete system of equations is given by:

$$\begin{aligned} {\left\{ \begin{array}{ll} \sum _{e=1}^{N_e}\left( {\mathbb {M}}_e\cdot \frac{{\mathbf {U}}^{k+1}_e-{\mathbf {U}}^{k}_e}{\tau }+{\mathbb {K}}_e\cdot {\mathbf {U}}^{k+1}_e-\gamma {\mathbf {F}}^{k+1}_e \right) ={\mathbf {0}},\\ \\ \sum _{e=1}^{N_e}\left( {\mathbb {M}}_e\cdot \frac{{\mathbf {V}}^{k+1}_e-{\mathbf {V}}^{k}_e}{\tau }+d{\mathbb {K}}_e\cdot {\mathbf {V}}^{k+1}_e-\gamma {\mathbf {G}}^{k+1}_e \right) ={\mathbf {0}}. \end{array}\right. } \end{aligned}$$

(A.9)

Notice that since the nonlinear terms ${\mathbf {F}}_e$ and ${\mathbf {G}}_e$ implicitly depend on time they are computed at $t^{k+1}$.

1.4 Appendix A.4: Fully Discrete Formulation

Let us now describe the numerical scheme to solve the nonlinear system of equations in Eq. (A.9). Let us rewrite Eq. (A.9) as:

$$\begin{aligned} {\left\{ \begin{array}{ll} \sum _{e=1}^{N_e}\left( {\mathbb {M}}_e\cdot \frac{{\mathbf {U}}^{k+1}_e-{\mathbf {U}}^{k}_e}{\tau }+{\mathbb {K}}_e\cdot {\mathbf {U}}^{k+1}_e-\gamma {\mathbf {F}}^{k+1}_e \right) ={\mathbf {R}}_U,\\ \\ \sum _{e=1}^{N_e}\left( {\mathbb {M}}_e\cdot \frac{{\mathbf {V}}^{k+1}_e-{\mathbf {V}}^{k}_e}{\tau }+d{\mathbb {K}}_e\cdot {\mathbf {V}}^{k+1}_e-\gamma {\mathbf {G}}^{k+1}_e \right) ={\mathbf {R}}_V, \end{array}\right. } \end{aligned}$$

(A.10)

where ${\mathbf {R}}_U$ and ${\mathbf {R}}_V$ are the residuals of the approximation. The Newton–Raphson scheme is then given by:

$$\begin{aligned} \begin{bmatrix} \dfrac{\partial {\mathbf {R}}_U}{\partial {\mathbf {U}}^{k+1}}&{}\dfrac{\partial {\mathbf {R}}_U}{\partial {\mathbf {V}}^{k+1}}\\ \dfrac{\partial {\mathbf {R}}_V}{\partial {\mathbf {U}}^{k+1}}&{}\dfrac{\partial {\mathbf {R}}_V}{\partial {\mathbf {V}}^{k+1}} \end{bmatrix} \begin{bmatrix} \Delta {\mathbf {U}}^{k+1}\\ \Delta {\mathbf {V}}^{k+1} \end{bmatrix} =- \begin{bmatrix} {\mathbf {R}}_U\\ {\mathbf {R}}_V \end{bmatrix}, \end{aligned}$$

(A.11)

where $\Delta {\mathbf {U}}^{k+1}$ and $\Delta {\mathbf {V}}^{k+1}$ are the differences of the approximations of ${\mathbf {U}}^{k+1}$ and ${\mathbf {V}}^{k+1}$, respectively, between solver iterations. The partial derivatives in Eq. (A.11) are given by:

$$\begin{aligned} \begin{aligned}&\dfrac{\partial {\mathbf {R}}_U}{\partial {\mathbf {U}}^{k+1}}=\sum _{e=1}^{N_e}\left( \frac{1}{\tau }{\mathbb {M}}_e+{\mathbb {K}}_e-\gamma \int _{\Omega ^h_e}\chi _i\chi _j \frac{\partial f}{\partial u^h}^{k+1} d{\Omega ^h_e} \right) ,\\&\dfrac{\partial {\mathbf {R}}_U}{\partial {\mathbf {V}}^{k+1}}=\sum _{e=1}^{N_e}\left( -\gamma \int _{\Omega ^h_e}\chi _i\chi _j \frac{\partial f}{\partial v^h}^{k+1} d{\Omega ^h_e} \right) ,\\&\dfrac{\partial {\mathbf {R}}_V}{\partial {\mathbf {U}}^{k+1}}=\sum _{e=1}^{N_e}\left( -\gamma \int _{\Omega ^h_e}\chi _i\chi _j \frac{\partial g}{\partial u^h}^{k+1} d{\Omega ^h_e} \right) ,\\&\dfrac{\partial {\mathbf {R}}_V}{\partial {\mathbf {V}}^{k+1}}=\sum _{e=1}^{N_e}\left( \frac{1}{\tau }{\mathbb {M}}_e+d{\mathbb {K}}_e-\gamma \int _{\Omega ^h_e}\chi _i\chi _j \frac{\partial g}{\partial v^h}^{k+1} d{\Omega ^h_e} \right) .\\ \end{aligned} \end{aligned}$$

(A.12)

Appendix B: Further Details, Figures and Tables

1.1 Appendix B.1: Wave Number Tables

See Tables 1 and 2.

Table 1 Combinations of d and $\gamma$ for the Schnakenberg reaction model with $a=0.1$ and $b=0.9$ for some one-dimensional wave modes

Table 2 Combinations of d and $\gamma$ for the Schnakenberg reaction model with $a=0.1$ and $b=0.9$ for some two-dimensional wave modes (Garzón-Alvarado 2007)

1.2 Appendix B.2: Changes in Polarity and Orientation

See Fig. 14.

Fig. 14

Schematic diagrams of pattern polarity and orientation differences. a A reference pattern with a change of polarity in b and a change in orientation in c

1.3 Appendix B.3: Schematic Representation of $c(k_2)$

See Fig. 15.

Fig. 15

Schematic representation of $c(k_2)$

1.4 Appendix B.4: Pattern Normalisation

$$\begin{aligned} v_i^*=\frac{v_i-\min (v)}{\max (v)-\min (v)}, \end{aligned}$$

(B.1)

where $v_i$ is the nodal value of v at the ith node (given by the finite element method discretisation).

1.5 Appendix B.5: Discretisation and Convergence Criteria

Regarding the discretisation, one-dimensional simulations were performed with uniform meshes of 1001 nodes and 1000 elements, and the two-dimensional simulations with uniform meshes of 10,000 elements and 10,201 nodes. Additionally, a time step of 0.01 was defined.

The discrete $L^2$-norm time derivative was calculated globally, as presented in Sarfaraz and Madzvamuse (2017):

$$\begin{aligned} \frac{||\mathbf {u}^{\mathbf {t}+{\varvec{\Delta }} \mathbf {t}}-\mathbf {u}^{\mathbf {t}||}}{\Delta t},\quad \text {and} \quad \frac{||\mathbf {v}^{\mathbf {t}+{\varvec{\Delta }} \mathbf {t}}-\mathbf {v}^{\mathbf {t}||}}{\Delta t}, \end{aligned}$$

(B.2)

and the simulations were stopped after a given tolerance, say $\varepsilon =10^{-6}$, was reached for both species.

1.6 Appendix B.6: Standard Convergence Graphics

See Fig. 16.

Fig. 16

Convergence for the unit one-dimensional domain. The bold line ( ) and the dashed line ( ) show the solutions of u and v, respectively. a Convergence to the heterogeneous solution with $(n_{\text {in}}=2,n_{\text {out}}=1)$ and $r=0.25$ and b convergence to the homogeneous solution with $(n_{\text {in}}=2,n_{\text {out}}=1)$ and $r=0.4$