A space-fractional-reaction-diffusion model for pattern formation in coral reefs

Abstract

In this paper, we propose a Space Fractional Reaction-Diffusion model for growth of corals in a tank. We analyse spatial temporal pattern formation behavior of the model through Turing type instability analysis. We determine the parameter regions of the model, in which the Turing instability occurs (the Turing instability space). We then discuss the effects of the fractional order and the model parameters on the spatial temporal patterns of the model. We investigate numerical solutions of the model using the Fourier Spectral method, two Euler type schemes and the MATLAB functionODE15s. The main advantage of using the Fourier spectral method is ability to represent the space fractional operator in fully diagonal form and ability to extend straightforwardly to two and three spatial dimensions. To find the numerical solutions with the other three methods, we transform the model equations into a system of ODEs by applying Matrix Transfer Technique and solve those system of ODEs. We compare the numerical solutions obtained by these methods and efficiencies of these methods.

Keywords:

• fractional-reaction-diffusion equations
• instability analysis
• patterns formations
• spectral methods

Public Interest Statement

Modelling coral morphogenesis processes has become a significant area of research. In this research paper we propose a mathematical model for pattern formation in coral reefs, applying physical phenomena called fractional diffusion. The phenomena fractional diffusion is applied instead of standard diffusion, as it is more suitable to model complex processes like pattern formation in coral reefs. We analytically proved that solutions of the proposed model form spatial temporal patterns when model parameters satisfy specific conditions. Under the same conditions, we solved the model and observed that these numerical solutions form branching structures like branching corals. We discussed the effects of the model parameters and the fractional exponent on the pattern formation behaviour.

1. Introduction

The purpose of this investigation is to explore the suitability of space fractional diffusion phenomena in modelling coral reefs pattern formations. Fractional diffusion equations arise when modelling anomalous (deviations from normal) diffusions features that are observed in many systems (Li & Wang, 2003; Metzler & Klafter, 2000; Santoro, de Paula, Lenzi, & Evangelista, 2011; Sun, Chen, & Chen, 2009; Upadhyaya, Rieu, Glazier, & Sawada, 2001). This introductory section provides brief overviews of corals, coral modelling and fractional reaction diffusion models.

Coral reefs are formed by calcium carbonate skeletons. Reef skeletons are built by tiny coral animals called polyps by depositing dissolved calcium carbonate ions on the existing reef surface that make up large coral colonies. The morphogenesis of coral is complex and that complexity define the beauty of the coral world. Various modeling approaches on coral morphogenesis processes have been proposed (Kaandorp, Blom et al., 2008; Kaandorp, Lowe, Frenkel, & Sloot, 1996; Kaandorp, Sloot et al., 2005; Merks, 2003a,2003b; Merks, Hoekstra, Kaandorp, & Sloot, 2003a,2003b; Mistr & Bercovici, 2003). A reaction diffusion type mathematical model for growth of corals in a tank has been proposed by Somathilake and Janak (2012,2014) based on the model that considers the nutrient polyps interaction (Mistr & Bercovici, 2003). In this paper, we propose a new Fractional-Reaction-Diffusion model that reduces to the Reaction-Diffusion model presented in Somathilake and Wedagedera (2014). Then we investigate the spatial temporal pattern formation behaviour of the proposed model.

Fractional calculus originated with the idea of a fractional derivative first raised by G. F. A. De L’Hôpital (1661–1704). However, in the past few decades Mathematicians have paid attention on fractional calculus and considerable developments have occured. In recent years there has been an increasing interest in fractional reaction diffusion equations in diverse and widespread fields of science and engineering (Benchohra, Henderson, Ntouyas, & Ouahab, 2008; Diethelm & Freed, 1999; Hilfer, 2000; Liu, Yang, & Turner, 2011; Sabatier, Agrawal, & Machado, 2007; Schiessel, Metzler, Blumen, & Nonnenmacher, 1995). Fractional reaction diffusion equations arise when modelling processes with non-locality and spatial heterogeneity. Especially, space fractional reaction diffusion models arise when modelling reaction and diffusion processes in heterogeneous media. Some analytical (Henry, Langlands, & Wearne, 2005; Jafari et al., 2013; Kirane, Ahmad, Alsaedi, & Al-Yami, 2014) and numerical (Bhrawy, Baleanu, & Mallawi, 2015; Bueno-Orovio, Kay, & Burrage, 2014; Bueno-Orovio & Pérez-Garca, 2006; Burrage, Hale, & Kay, 2012; Ilić, Liu, Turner, & Anh, 2005,2006; Pindza & Owolabi, 2016; Rida, El-Sayed, & Arafa, 2010; Tchier, Inc, Korpinar, & Baleanu, 2016; Yang, Liu, & Turner, 2010; Yu, Liu, Anh, & Turner, 2008; Zhuang, Liu, Anh, & Turner, 2009) techniques play an important role in identifying the solution behaviour and obtaining approximate solutions of fractional differential equations.

The general form of the one variable fractional reaction diffusion equation is:(1) $$\begin{array}{c}\hfill {\displaystyle \frac{\partial u}{\partial t}}=-D_{\gamma }{(-\mathrm{\Delta })}^{\gamma /2}u+f(u,t).\end{array}$$(1)

Here $-D_{\gamma }{(-\mathrm{\Delta })}^{\gamma /2}u$, $D_{\gamma }$ and $f(u,v)$ represent a fractional diffusion process, diffusion rate and the reaction process, respectively.

For $\gamma =2$, anomalous diffusion reduces to normal diffusion. For $\gamma <2$$(\gamma >2)$, the diffusion process is slower (faster) than normal diffusion, where it is called sub-diffusive (super-diffusive). Sub-diffusion terms arise when we model processes related to spatial or temporal constraints such as fractured and porous media and fractal lattices (Henry & Wearne, 2000). Super-diffusion terms may arise when we model processes such as particles transportation through chaotic or turbulent media (Henry & Wearne, 2000). Coral reef mass can be considered as a porous medium. Therefore, in order to make the model more realistic, anomalous diffusion is considered in the model formation process.

In order to allow spatial temporal patterns to form by the solutions of a reaction diffusion system, the solutions should be stable in the absence of the diffusion, and be unstable in the presence of the diffusion (Murray, 2003; Turing, 1952). Similar arguments have been used in space fractional reaction diffusions equations (Henry et al., 2005; Tian, Zhou, & Deng, 2015). In this paper we perform a Turing type instability analysis of the proposed fractional reaction diffusion model and we determine instability regions (parameter ranges in which Turing instability occurs). Then we investigate effects of the fractional exponent on the pattern formation behaviour.

This paper is organized as follows. An explanation on the derivation of the mathematical model is given in Section 2. Turing type instability analysis of the model is performed in Section 3. Also, in the same section, we discuss the effects of the fractional exponent of the spatial diffusion and model parameters on spatial pattern formation by the solutions of the model. Some numerical methods based on finite difference and Fourier spectral methods are explained in Section 4. In Section 5, numerical results on 1D and 2D spaces are presented. Finally, we present the discussion and conclusions based on the results in Section 6.

2. Mathematical model

In this section we derive a Fractional-Reaction-Diffusion model for growth of coral in a tank. Consider a water filled tank with some coral larva (planula) settled on the bottom of the tank. Assume that the coral growth is taking place in this tank, and there is no fluid motion within the tank. Also, assume that coral polyps consume dissolved nutrients and produce additional solid material (corals produce their skeleton by extracting dissolved calcium carbonate surrounding them). This process can be regarded as a reaction process between dissolved nutrient and dissolved calcium carbonate (Mistr & Bercovici, 2003). Let A and B denote the dissolved nutrients and dissolved solid material (calcium carbonate), respectively, and u and v denote their respective biomasses. In this case, nutrient and solid material concentrations at a point is controlled by two processes, reaction and anomalous diffusion.

Thus, nutrients react with the dissolved solid material (calcium carbonate ions) and produce additional solid material. Assume that the reaction process of A and B takes the form (as in Mistr & Bercovici, 2003):$$A+2B\stackrel{b_1}{⟶}3B.$$

Assume that nutrients are supplied to the system in a constant rate and portion of the nutrients waste proportional to its instantaneous concentration. Then we can describe the time evolution process of nutrients as follows:$$\begin{array}{cc}& \left(\text{Time rate change of Nutrient concentration in reaction process}\right)=\left(\text{Anomalous diffusion of}u\right)\hfill \\ \hfill & +\left(\text{Constant supply rate of}u\right)-\left(\text{wasting rate of}u\right)-\left(\text{Reactive loss of}u\right).\hfill \end{array}$$

Assume that constant rates of nutrients supplying, wasting and reaction are a, b, and $b_1$ respectively. Let $d_u$ and $d_v$ be the diffusion rates of dissolved nutrient and calcium carbonate respectively. Then we can directly write the rate equation for the above process as:(2) $$\begin{array}{c}\hfill \frac{\partial u}{\partial t}=-d_u{(-\mathrm{\Delta })}^{\gamma /2}u+a-bu-b_{1}u{v}^2;1<\gamma \le 2.\end{array}$$(2)

This system can be written in the form:(3) $$\begin{array}{c}\hfill {\displaystyle \frac{\partial u}{\partial t}}=-d_u{(-\mathrm{\Delta })}^{\gamma /2}u+b(u_s-u)-b_1{v}^{2}u\end{array}$$(3)

where $u_s={\displaystyle \frac{a}{b}}$. That is $b(u-u_s)$ is the overall nutrient supply rate to the system.

It is reasonable to assume that the time rate change of the dissolved solid material concentration is controlled by fractional diffusion, reactive production and loss due to aggregation on existing coral reefs. This process can be represented in the form:$$\begin{array}{cc}& \left(\text{Time rate change of dissolved calcium carbonate ions concentration}v\right)=\left(\text{Anomalous diffusion of}v\right)\hfill \\ \hfill & -\left(\text{Loss of}v\text{due to deposition}\right)+\left(\text{Reactive production of}v\right)\hfill \end{array}$$

which gives(4) $$\begin{array}{c}\hfill {\displaystyle \frac{\partial v}{\partial t}}=-d_v{(-\mathrm{\Delta })}^{\gamma /2}v-b_{2}v+b_1{v}^{2}u,\end{array}$$(4)

where $b_2$ is the depositing rate of calcium carbonate on the existing corals.

Let w denotes the depositing calcium carbonate concentration on the coral reef. This process can be represented by the equation:(5) $$\begin{array}{c}\hfill {\displaystyle \frac{\partial w}{\partial t}}=b_{2}v.\end{array}$$(5)

Finally we get a set of fractional reaction diffusion equations(6) $$\begin{array}{cc}\hfill {\displaystyle \frac{\partial u}{\partial t}}& =-d_u{(-\mathrm{\Delta })}^{\gamma /2}u+b(u_s-u)-b_1{v}^{2}u\hfill \\ \hfill {\displaystyle \frac{\partial v}{\partial t}}& =-d_v{(-\mathrm{\Delta })}^{\gamma /2}v-b_{2}v+b_1{v}^{2}u\hfill \\ \hfill {\displaystyle \frac{\partial w}{\partial t}}& =b_{2}v.\hfill \end{array}$$(6)

There are six parameters in this model. In order to reduce the number of parameters, we use nondimensionalisation techniques.

2.1. Nondimensionalisation

Now we nondimensionalise the equations system (6) by applying the coordinate transformations $t^{\ast }=bt$, $x^{\ast }=\frac{x}{L}$, $y^{\ast }=\frac{y}{L}$, $z^{\ast }=\frac{z}{L}$, where $L={\left({\displaystyle \frac{d_u}{b}}\right)}^{1/\gamma }$, we get(7) $$\begin{array}{cc}\hfill {\displaystyle \frac{\partial \overline{u}}{\partial t^{\ast }}}& ={-(-{\mathrm{\Delta }}^{\ast })}^{\gamma /2}\overline{u}+1-\overline{u}-{\alpha }^2{\overline{v}}^2\overline{u}\hfill \\ \hfill {\displaystyle \frac{\partial \overline{v}}{\partial t^{\ast }}}& =-d{(-{\mathrm{\Delta }}^{\ast })}^{\gamma /2}\overline{v}-\lambda \overline{v}+{\alpha }^2{\overline{v}}^2\overline{u}\hfill \\ \hfill {\displaystyle \frac{\partial \overline{w}}{\partial t^{\ast }}}& =\lambda \overline{v},\hfill \end{array}$$(7)

where $\overline{u}={\displaystyle \frac{u}{u_s}}$, $\overline{v}={\displaystyle \frac{v}{u_s}}$, ${\alpha }^2={\displaystyle \frac{b_1{u}_s^2}{b}}$ and $\lambda ={\displaystyle \frac{b_2}{b}}$, $d={\displaystyle \frac{d_v}{d_u}}$.

For notational convenience, eliminating bars and stars we obtain(8) $$\begin{array}{cc}\hfill {\displaystyle \frac{\partial u}{\partial t}}& =-{(-\mathrm{\Delta })}^{\gamma /2}u+1-u-{\alpha }^{2}u{v}^2\hfill \\ \hfill {\displaystyle \frac{\partial v}{\partial t}}& =-d{(-\mathrm{\Delta })}^{\gamma /2}v-\lambda v+{\alpha }^{2}u{v}^2\hfill \\ \hfill {\displaystyle \frac{\partial w}{\partial t}}& =\lambda v.\hfill \end{array}$$(8)

2.2. Boundary conditions and Matrix Transfer technique

As we consider the coral growth in a tank, homogeneous Neumann boundary conditions (also referred to as zero flux or insulating boundary conditions) are more reasonable for standard diffusion ($\gamma =2$). That is, the boundary conditions can be expressed as:(9) $$\begin{array}{c}\hfill \left.\begin{array}{cc}\mathrm{\nabla }u·\mathbf{n}=0,\hfill & x\in \partial \mathrm{\Omega }\hfill \\ \mathrm{\nabla }v·\mathbf{n}=0,\hfill & x\in \partial \mathrm{\Omega }\hfill \end{array}\right\},\end{array}$$(9)

where $\mathbf{n}$ denotes the outwards unit normal vector to the boundary $\partial \mathrm{\Omega }$.

However, the fractional diffusion operator representing anomalous diffusion is a non local operator, and implementation of the problems on finite domains with boundary conditions except zero Dirichlet boundary conditions is a challenge. For the case of zero Dirichlet boundary conditions, fractional differential equation can be extended, so that its solution is equal to zero on the complement of the bounded domain (Yang et al., 2010).

In handling fractional reaction diffusion problems on bounded domains, the method called Matrix Transfer Technique, proposed by Ilić et al. (2005,2006) plays an important role.

Let A be the discrete Laplacian coupled with standard homogeneous Neumann boundary conditions and let $A=VD{V}^{-1}$, where V is the matrix, of which columns are eigenvectors of A, and D is the diagonal matrix of which elements are eigenvalues of A. The discrete Laplacian operator, A, on $[0,L]$ coupled with homogeneous Neumann boundary conditions, obtained with a finite difference approximation on a uniform mesh of $n+1$ nodes with step size $h=L/n$ and with a second order approximation of the boundary conditions is given by:(10) $$\begin{array}{c}\hfill A={\frac{1}{h^2}\left(\begin{array}{cccc}2\hfill & -2\hfill & \hfill & \\ -1\hfill & 2\hfill & -1\hfill & \\ \hfill & \ddots \hfill & \ddots \hfill & \ddots \hfill & \\ \hfill & \hfill & -1\hfill & 2\hfill & -1\hfill \\ \hfill & \hfill & -2\hfill & 2\hfill \end{array}\right).}_{(n+1)\times (n+1)}\end{array}$$(10)

According to Ilić et al. (2005,2006), the discrete fractional operator $A^{\gamma /2}$ ($1<\gamma \le 2$) is defined as $A^{\gamma /2}=V{D}^{\gamma /2}{V}^{-1}$. This formulation is called the Matrix Transfer Technique. There is an issue whether the boundary conditions embedded in A transfer properly to the fractional order of the matrix, $A^{\gamma /2}$. Very recently, it has been proved that Neumann boundary conditions embedded in A transfers properly to $A^{\gamma /2}$ in the one dimensional case (Cusimano, Burrage, Turner, & Kay, 2017).

3. Turing type instability of the model

In this section, considering the linear stability of the steady state, we present the Turing parameter space of the model. We derive the conditions for Turing bifurcation through linear stability analysis of the system (8).

3.1. Steady states

Now we represent the first two equations of the system (8) of the form:(11) $$\begin{array}{cc}\hfill {\displaystyle \frac{\partial u}{\partial t}}& =-{(-\mathrm{\Delta })}^{\gamma /2}u+F(u,v),\hfill \\ \hfill {\displaystyle \frac{\partial v}{\partial t}}& =-d{(-\mathrm{\Delta })}^{\gamma /2}v+G(u,v),\hfill \end{array}$$(11)

where $F(u,v)=1-u-{\alpha }^{2}u{v}^2$ and $G(u,v)=-\lambda v+{\alpha }^{2}u{v}^2$. There are three homogeneous steady states $S_1\equiv (u_{s1},v_{s1})$, $S_2\equiv (u_{s2},v_{s2})$ and $S_3\equiv (u_{s3},v_{s3})$ for the corresponding system of ordinary differential equations. Here $u_{s1}=1$, $v_{s1}=0$, $u_{s2}={\displaystyle \frac{\alpha -\sqrt{{\alpha }^2-4{\lambda }^2}}{2\alpha }}$, $v_{s2}={\displaystyle \frac{\alpha +\sqrt{{\alpha }^2-4{\lambda }^2}}{2\alpha \lambda }}$, $u_{s3}={\displaystyle \frac{\alpha +\sqrt{{\alpha }^2-4{\lambda }^2}}{2\alpha }}$ and $v_{s3}={\displaystyle \frac{\alpha -\sqrt{{\alpha }^2-4{\lambda }^2}}{2\alpha \lambda }}$ for $\alpha >2\lambda$.

The Jacobian matrix evaluated at $S_i$ is(12) $$\begin{array}{c}\hfill A_i={\left.\left(\begin{array}{cc}{F}_u\hfill & F_v\hfill \\ G_u\hfill & G_v\hfill \end{array}\right)\right|}_{S_i}={\left.\left(\begin{array}{cc}{a}_{11}\hfill & a_{12}\hfill \\ a_{21}\hfill & a_{22}\hfill \end{array}\right)\right|}_{S_i}=\left(\begin{array}{cc}-1-{\alpha }^2{v}_{\mathit{si}}^2& -2\lambda \\ {\alpha }^2{v}_{\mathit{si}}^2& \lambda \end{array}\right).\end{array}$$(12)

The trace and determinant of A are $Tr(A_i)=\lambda -1-{\alpha }^2{v}_{\mathit{si}}^2$ and $Det(A_i)=\lambda ({\alpha }^2{v}_{\mathit{si}}^2-1)$, respectively.

3.2. Instability conditions

We can show that the trivial steady state, $S_1$ is a stable node, $S_3$ is a saddle point and the stability of $S_2$ further depends on the parameters $\alpha$ and $\lambda$. We note the above conditions are not satisfied at $S_1$. Since $S_3$ is a saddle point, Turing instability does not hold at $S_3$. Therefore, a Turing type instability occurs only at $S_2$ and we focus on the Turing type patterns at $S_2$.

Some reaction diffusion systems have special characteristics. They have stable steady states in the absence of diffusion (in a spatially homogeneous system) and these states become unstable in the presence of diffusion (in a spatially inhomogeneous system). Solutions of such systems form spatial temporal patterns. Alan Turing explained this phenomena for reaction diffusion models arising in chemistry (Turing, 1952) and these conditions are called Turing’s instability conditions. A similar explanation can be done for fractional reaction diffusion equations. That is, fractional reaction diffusion systems may form spatial patterns when a Turing type instability occurs. J. D. Murray has obtained Turing instability conditions for the following general form of two components reaction diffusion system subject to homogeneous zero flux boundary conditions (Murray, 2003).$$\begin{array}{c}\hfill \left.\begin{array}{cc}{\displaystyle \frac{\partial u}{\partial t}}=\mathrm{\Delta }u+f(u,v),x\in \mathrm{\Omega }\subset {\mathbb{R}}^2,\hfill & t>0\hfill \\ {\displaystyle \frac{\partial v}{\partial t}}=d\mathrm{\Delta }v+g(u,v),x\in \mathrm{\Omega }\subset {\mathbb{R}}^2,\hfill & t>0.\hfill \end{array}\right\}\end{array}$$

Here, $f(u,v)$, and $g(u,v)$ are functions of u and v, representing the reaction processes of the relevant model, and $\mathrm{\Omega }$ is the domain of the considered model. The Turing instability conditions of this model take the form:$$\begin{array}{c}\hfill C_1\equiv \text{Det}(A)>0,\\ \hfill C_2\equiv \text{Tr}(A)<0,\\ \hfill C_3\equiv d{a}_{11}+a_{22}>0,\\ \hfill C_4\equiv {C_3}^2-4d{C}_1>0;\\ \hfill \end{array}$$

where$$\begin{array}{c}\hfill A={\left.\left(\begin{array}{cc}{f}_u\hfill & f_v\hfill \\ g_u\hfill & g_v\hfill \end{array}\right)\right|}_S={\left.\left(\begin{array}{cc}{a}_{11}\hfill & a_{12}\hfill \\ a_{21}\hfill & a_{22}\hfill \end{array}\right)\right|}_S\end{array}$$

and $s\equiv (u_0,v_0)$ is a stable steady state of the corresponding spatially homogeneous system. Now we derive Turing type instability conditions for the system (11). By linearising the system (11) about a steady state $(u_0,v_0)$ we get$$\begin{array}{c}\hfill \left.\begin{array}{c}{\displaystyle \frac{\partial u_1}{\partial t}}=-{(-\mathrm{\Delta })}^{\gamma /2}{u}_1+F_u(u_0,v_0)u_1+F_v(u_0,v_0)v_1\hfill \\ {\displaystyle \frac{\partial v_1}{\partial t}}=-d{(-\mathrm{\Delta })}^{\gamma /2}{v}_1+G_u(u_0,v_0)u_1+G_v(u_0,v_0)v_1\hfill \end{array}\right\}\end{array}$$

where $u_1$ and $v_1$ are small perturbations of $(u_0,v_0)$ such that $u=u_0+u_1$ and $v=v_0+v_1$. The matrix form of this system is:(13) $$\begin{array}{c}\hfill \left(\begin{array}{c}{\displaystyle \frac{\partial u_1}{\partial t}}\hfill \\ {\displaystyle \frac{\partial v_1}{\partial t}}\hfill \end{array}\right)=\left(\begin{array}{cc}-{(-\mathrm{\Delta })}^{\gamma /2}+a_{11}& a_{12}\\ a_{21}& a_{22}-d{(-\mathrm{\Delta })}^{\gamma /2}\end{array}\right)\left(\begin{array}{c}{u}_1\hfill \\ v_1\hfill \end{array}\right),\end{array}$$(13)

where $a_{11}=F_u(u_0,v_0)$, $a_{12}=F_v(u_0,v_0)$, $a_{21}=G_u(u_0,v_0)$ and $a_{22}=G_v(u_0,v_0)$. Consider the solutions of the form:(14) $$\begin{array}{c}\hfill \left(\begin{array}{c}{u}_1\hfill \\ v_1\hfill \end{array}\right)=\left(\begin{array}{c}A\hfill \\ B\hfill \end{array}\right)e^{\sigma t+i\mathbf{k}·\mathbf{r}}.\end{array}$$(14)

Substituting this solution in equation (13) we obtain:(15) $$\begin{array}{c}\hfill \left(\begin{array}{c}\sigma u_1\hfill \\ \sigma v_1\hfill \end{array}\right)=\left(\begin{array}{cc}{a}_{11}-k^{\gamma }& a_{12}\\ a_{21}& a_{22}-d{k}^{\gamma }\end{array}\right)\left(\begin{array}{c}{u}_1\hfill \\ v_1\hfill \end{array}\right).\end{array}$$(15)

From this, we get the dispersion relation(16) $$\begin{array}{c}\hfill \left|\begin{array}{cc}\sigma -a_{11}+k^{\gamma }\hfill & -a_{12}\hfill \\ -a_{21}\hfill & \sigma +d{k}^{\gamma }-a_{22}\hfill \end{array}\right|=0,\end{array}$$(16)

where $k=\Vert \mathbf{k}\Vert$. This can be written in the form(17) $$\begin{array}{c}\hfill {\sigma }^2-g(\mu )\sigma +h(\mu )=0.\end{array}$$(17)

Here $\mu =k^{\gamma }$, $g(\mu )=a_{11}+a_{22}-\mu -d\mu$ and $h(\mu )=d{\mu }^2-(d{a}_{11}+a_{22})\mu +a_{11}{a}_{22}-a_{12}{a}_{21}$. In terms of the model parameters$$\begin{array}{cc}\hfill g(\mu )& =\frac{\left({\alpha }^2-2{\lambda }^3+\alpha \sqrt{{\alpha }^2-4{\lambda }^2}+2{\lambda }^2\mu +2d{\lambda }^2\mu \right)}{2{\lambda }^2}\text{and}\hfill \\ \hfill h(\mu )& =\frac{{\alpha }^2\left(\lambda +d\mu \right)+\alpha \sqrt{{\alpha }^2-4{\lambda }^2}\left(\lambda +d\mu \right)-2{\lambda }^2\left(-d{\mu }^2+\lambda \left(2+\mu \right)\right)}{2{\lambda }^2}.\hfill \end{array}$$

Since $S_2$ is a stable steady state, $C_1^{\prime }\equiv a_{11}{a}_{22}-a_{12}{a}_{21}>0$, $C_2^{\prime }\equiv a_{11}+a_{22}<0$ and hence $g(\mu )<0$. Therefore, the real parts of the solutions for $\sigma$ in equation (13) are positive only if $h(\mu )<0$. That is, the system undergoes Turing bifurcation if and only if $h(\mu )<0$. Now we find the critical values of $h(\mu )$.$$\begin{array}{cc}& {\displaystyle \frac{dh(\mu )}{d\mu }}=2d\mu -(d{a}_{11}+a_{22})\hfill \\ \hfill & {\displaystyle \frac{dh(\mu )}{d\mu }}=0\Rightarrow 2d\mu -(d{a}_{11}+a_{22})=0\Rightarrow \mu ={\mu }_c=k_c^{\gamma }={\displaystyle \frac{d{a}_{11}+a_{22}}{2d}}.\hfill \end{array}$$

Then $h(\mu )$ has a single minimum, $h_{\text{min}}(\mu )={\displaystyle \frac{4d(a_{11}{a}_{22}-a_{12}{a}_{21})-{(d{a}_{11}+a_{22})}^2}{4d}}$, at $\mu ={\mu }_c=k_c^{\gamma }$. In order for spatial patterns to exist, $h_{\text{min}}<0$ and $k_c>0$. That is the conditions: $C_3^{\prime }\equiv d{a}_{11}+a_{22}>0$ and $C_4^{\prime }\equiv 4d(Det(A))-{(d{a}_{11}+a_{22})}^2<0$

should be satisfied. It is clear that $C_i=C_i^{\prime }$, for $i=1,2,3,4.$

Therefore, the instability conditions for the fractional reaction diffusion equations, and for the standard Reaction Diffusion equations are same. But the growth rate, $\sigma$, and wave numbers, k, depend on $\gamma$. In terms of the model parameters, we can represent the instability conditions at $S_2$ in the following form:$$\begin{array}{cc}& C_1\equiv {\displaystyle \frac{{\alpha }^2-4{\lambda }^2+\alpha \sqrt{{\alpha }^2-4{\lambda }^2}}{2\lambda }}>0,\hfill \\ \hfill & C_2\equiv -{\displaystyle \frac{{\alpha }^2-2{\lambda }^3+\alpha \sqrt{{\alpha }^2-4{\lambda }^2}}{2{\lambda }^2}}<0,\hfill \\ \hfill & C_3\equiv {\displaystyle \frac{2{\lambda }^3-\alpha d(\alpha +\sqrt{{\alpha }^2-4{\lambda }^2})}{2{\lambda }^2}}>0,\hfill \\ \hfill & C_4\equiv {\left({\displaystyle \frac{2{\lambda }^3-d\alpha (\alpha +\sqrt{{\alpha }^2-4{\lambda }^2})}{2{\lambda }^2}}\right)}^2-4d\left({\displaystyle \frac{{\alpha }^2-4{\lambda }^2+\alpha \sqrt{{\alpha }^2-4{\lambda }^2}}{2\lambda }}\right)>0.\hfill \end{array}$$

3.3. Instability regions

It is clear that $C_1$ holds at $S_2$ in the region $\alpha >2\lambda$. The bifurcation $C_2=0$ (Hopf Bifurcation), gives one positive real roots for $\alpha$ when $\lambda >1$. These are$${\alpha }_1^{(2)}={\displaystyle \frac{{\lambda }^2}{\sqrt{\lambda -1}}}.$$

Equation $C_3=0$ gives two real roots for $\alpha$ when $\lambda >d$. These are$${\alpha }_{1,2}^{(3)}={\displaystyle \frac{\pm {\lambda }^2}{\sqrt{d(\lambda -d)}}}.$$

Equation $C_4=0$ gives two positive real roots for $\alpha$ when $\lambda >d$. These are$${\alpha }_{1,2}^{(4)}={\displaystyle \frac{\sqrt{d{\lambda }^3(8d^2+7d\lambda +3{\lambda }^2)\mp 2\sqrt{2}d\sqrt{{\lambda }^7(\lambda -d){(2d+\lambda )}^2}}}{d(d+\lambda )}}.$$

Since $\alpha >2\lambda$, letting $\alpha =2\lambda +ϵ$ ($ϵ>0$) in $C_2$, we obtain(18) $$\begin{array}{cc}\hfill C_2& =-{\displaystyle \frac{{(2\lambda +ϵ)}^2-2{\lambda }^3+(\text{positive term})}{2{\lambda }^2}}\hfill \\ \hfill & =-{\displaystyle \frac{2{\lambda }^2(2-\lambda )+(\text{positive term})}{2{\lambda }^2}}\hfill \\ \hfill & <0;\text{if}\lambda \le 2.\hfill \end{array}$$(18)

Therefore $C_2$ is negative for $\lambda \le 2$. Also, ${\alpha }_1^{(3)}>{\alpha }_1^{(4)}$ for $\lambda \ne 2d$ and ${\alpha }_1^{(3)}={\alpha }_1^{(4)}$ for $\lambda =2d$. Moreover(19) $$\begin{array}{cc}\hfill C_3& =-{\displaystyle \frac{d{(2\lambda +ϵ)}^2-2{\lambda }^3+(\text{positive term})}{2{\lambda }^2}}\hfill \\ \hfill & =-{\displaystyle \frac{2{\lambda }^2(2d-\lambda )+(\text{positive term})}{2{\lambda }^2}}\hfill \\ \hfill & <0;\text{if}\lambda \le 2d.\hfill \end{array}$$(19)

Therefore, $C_3$ is not satisfied in the region $\lambda \le 2d$. Also, we can show numerically that

• $2\lambda <{\alpha }_1^{(4)}<{\alpha }_1^{(3)}<{\alpha }_2^{(4)}$ for $2d<\lambda <2$ and ${\alpha }_1^{(3)}={\alpha }_1^{(4)}$ at $\lambda =2d$.

• $2\lambda <{\alpha }_1^{(2)}<{\alpha }_1^{(4)}<{\alpha }_1^{(3)}<{\alpha }_2^{(4)}$ for $2<\lambda$ and ${\alpha }_1^{(2)}={\alpha }_1^{(4)}$ at $\lambda =2$.

• $C_1$ is satisfied in the region $\{\alpha \in \mathbb{R}|2\lambda <\alpha \}$.

• $C_2$ is satisfied in the region $\{\alpha \in \mathbb{R}|\alpha >{\alpha }_1^{(2)},\lambda >2\}\cup \{\alpha \in \mathbb{R}|\lambda <2\}$.

• $C_3$ is satisfied in the region $\{\alpha \in \mathbb{R}|\alpha <{\alpha }_1^{(3)},2d<\lambda \}$.

• $C_4$ is satisfied in the region $\{\alpha \in \mathbb{R}|\alpha <{\alpha }_1^{(4)},d<\lambda \}\cup \{\alpha \in \mathbb{R}|\alpha >{\alpha }_2^{(4)},d<\lambda \}$.

By considering these facts, we can write the instability region (Turing space) as a union of two regions ${\mathcal{R}}_1$ and ${\mathcal{R}}_2$. Here,$$\begin{array}{cc}& {\mathcal{R}}_1=\{(\lambda ,\alpha );2d\le \lambda \le 2,2\lambda <\alpha <{\alpha }_{\text{max}}(\lambda ,d)\},\hfill \\ \hfill & {\mathcal{R}}_2=\{(\lambda ,\alpha );\lambda >2,{\alpha }_{\text{min}}(\lambda )<\alpha <{\alpha }_{\text{max}}(\lambda ,d)\},\hfill \end{array}$$

where ${\alpha }_{\text{min}}(\lambda )={\alpha }_1^{(2)}$ and ${\alpha }_{\text{max}}(\lambda ,d)={\alpha }_1^{(4)}$. Figure 1 shows the parameter regions ${\mathcal{R}}_1$ and ${\mathcal{R}}_2$ in $\alpha$, $\lambda$, d space as volumes.

Figure 1. Turing space for: (a) Region ${\mathcal{R}}_1$, (b) Region ${\mathcal{R}}_2$ (for $0<d<d_{c1}$), (c) Region ${\mathcal{R}}_2$ (for $d>d_{c1}$).

3.4. Critical values of d

When $d<d_{c_1}(\equiv 3\pm 2\sqrt{2})$ and $\lambda >2$, the instability parameter region ${\mathcal{R}}_2$ is unbounded and when $d>d_{c_1}$ and $\lambda >2$ it is bounded (Somathilake & Wedagedera, 2014). In order to satisfy Turing instability conditions the diffusion rate d should satisfy $d<d_c(\alpha ,\lambda )$ (Somathilake & Wedagedera, 2014). Here $d_c(\alpha ,\lambda )=\frac{3{\alpha }^2{\lambda }^3+3\alpha {\lambda }^3\sqrt{{\alpha }^2-4{\lambda }^2}-\mathcal{A}}{{\alpha }^2\left({\alpha }^2-2{\lambda }^2+\alpha \sqrt{{\alpha }^2-4{\lambda }^2}\right)}$,

where $\mathcal{A}=4\left(2{\lambda }^5+\sqrt{{\lambda }^6\left({\alpha }^4-5{\alpha }^2{\lambda }^2+4{\lambda }^4+{\alpha }^3\sqrt{{\alpha }^2-4{\lambda }^2}-3\alpha {\lambda }^2\sqrt{{\alpha }^2-4{\lambda }^2}\right)}\right).$

3.5. Effect of the fractional exponent and model parameters for the growth

Figures 2 and 3 show the variation of the real part of $\sigma$, Re($\sigma$), with respect to k at different levels of $\gamma$, and different levels of d, respectively for particular parameter values $\alpha =4.24$, $\lambda =2.1$ lying in region ${\mathcal{R}}_2$. As time increases, the solution of the system is dominated by the wave numbers corresponding to the maximum value of growth rate $\sigma$.

Figure 2. The variation of the real part of $\sigma$, Re($\sigma$), against k at different levels of $\gamma$ for particular d. Note: Other parameters which lie in region ${\mathcal{R}}_2$ are $\alpha =4.24$ and $\lambda =2.1$.

Consider the dispersion relation (17). Critical values of $\sigma$ are given by(20) $$\begin{array}{c}\hfill {\displaystyle \frac{d\sigma }{d\mu }}=0.\end{array}$$(20)

Solving (20) for $k^{\gamma }$ we get two solutions $k_1^{\gamma }$ and $k_2^{\gamma }$ given by$$\begin{array}{cc}& {\mu }_1=k_1^{\gamma }=\frac{A_2+B_2}{2d{\left(d-1\right)}^2{\lambda }^4},\hfill \\ \hfill & {\mu }_2=k_2^{\gamma }=\frac{A_2-B_2}{2d{\left(d-1\right)}^2{\lambda }^4},\hfill \end{array}$$

where$$\begin{array}{cc}\hfill A_2& =-d{\lambda }^2\left({\alpha }^2+2{\lambda }^3+\alpha \sqrt{{\alpha }^2-4{\lambda }^2}\right)+d^2{\lambda }^2\left({\alpha }^2+2{\lambda }^3+\alpha \sqrt{{\alpha }^2-4{\lambda }^2}\right)\text{and}\hfill \\ \hfill B_2& =2\sqrt{d{\left(-1+d^2\right)}^2{\lambda }^7\left({\alpha }^2-2{\lambda }^2+\alpha \sqrt{{\alpha }^2-4{\lambda }^2}\right).}\hfill \end{array}$$

We can show that $\sigma$ is maximized at $k=k_1$. That is ${\sigma }_{\text{max}}{(\alpha ,\lambda )=\sigma |}_{k=k_1}$. The plots of ${\sigma }_{\text{max}}$ against $\alpha$ when $\lambda =2.1$ for different levels of d are shown in Figure 4.

Figure 3. (a) The variation of Re($\sigma$) against k at different levels of d for particular $\gamma$ values. Note: The other parameters that lie in region ${\mathcal{R}}_2$ are $\alpha =4.24$, $\lambda =2.1$.

Figure 4. (a) The variation of ${\sigma }_{\text{max}}$ against $\alpha$ for different values of d when $\lambda =2.1$.

It is clear that ${\sigma }_{\text{max}}(\alpha ,\lambda )$ is maximized at $\alpha ={\alpha }_{\text{min}}(\lambda )$ for given $\lambda$. The wavelength corresponding to the maximum growth rate is given by ${\omega }_{\text{max}}=2\pi /k_1$. The wave number $n_{\text{max}}$ corresponding to the maximum growth rate is given by $n_{\text{max}}=k_{1}L/\pi$. The plot of $n_{\text{max}}$ against $\alpha$ for $\lambda =2.1$, $d=0.01$, at different levels of $\gamma$ is shown in Figure 5 . We can observe that for fixed $\lambda$ and d, the value of $n_{\text{max}}$ increases when $\gamma$ increases from 1. That is the heterogeneity of the spatial patterns increases as $\gamma$ increases when $\alpha$, $\lambda$ and d are fixed. Also we observe that when both $\lambda$ and $\alpha$ increase the possible minimum and maximum values of $n_{\text{max}}$ as well as the range of $n_{\text{max}}$ increase. This implies that the heterogeneity of the spatial patterns increases as $\alpha$ and $\lambda$ increase when d and $\gamma$ are fixed.

Figure 5. The variation of $n_{\text{max}}$ with respect to $\alpha$ for different levels of $\gamma$. Note: The values of other parameters are $\lambda =2.1$, $d=0.01$, respectively.

4. Numerical methods

4.1. Fourier spectral method

The Fourier spectral method has been used in numerical computations of fractional reaction diffusion equations by several researchers (Bueno-Orovio et al., 2014; Bueno-Orovio & Pérez-Garca, 2006; Pindza & Owolabi, 2016). The fractional Laplace operator on a bounded domain is expressed in Bueno-Orovio et al. (2014) and Ilić et al. (2006) based on a set of orthonormal eigenfunctions and corresponding eigenvalues as follows:

Definition 4.1

(Ilić et al., 2006) Suppose the Laplacian ($-\mathrm{\Delta }$) has a complete set of orthonormal eigenfunctions $\{{\varphi }_i\}$ corresponding eigenvalues ${\lambda }_i$ on a bounded domain $\mathcal{D}$, i.e. $(-\mathrm{\Delta }){\varphi }_i={\lambda }_i{\varphi }_i$ on $\mathcal{D}$; $\mathcal{B}(\varphi )=0$ on $\partial D$, where $\mathcal{B}(\varphi )$ is one of the three standard boundary conditions . Let$$U_{\beta }=\left\{u=\sum _{i=1}^{\infty }\widehat{u_i}{\varphi }_i,|\widehat{u_i}=<u,{\varphi }_i>,\sum _{i=1}^{\infty }|\widehat{u_i}{|}^2{|{\lambda }_i|}^{\beta }<\infty ,\beta =\text{max}(\gamma ,0)\right\}.$$

Then, for any $u\in U_{\beta }$, ${(-\mathrm{\Delta })}^{\gamma /2}$ is defined by$${(-\mathrm{\Delta })}^{\gamma /2}u=\sum _{i=1}^{\infty }\widehat{u_i}{{\lambda }_i}^{\gamma }{\varphi }_i,$$

here ${\lambda }_i$ and ${\varphi }_i$ depend on the specified boundary conditions.

The boundary conditions of our model corresponding to standard diffusion are homogeneous Neumann type, and according to Cusimano et al. (2017) these boundary conditions properly transfer to the fractional diffusion operator in one dimensional space. Therefore, the eigenvalues and eigenvectors of the fractional diffusion operator on a bounded domain in one dimensional space with homogeneous Neumann boundary conditions are well defined. In this paper we assume that the result holds for two dimension as well.

The eigenvalues ${\mu }_i$ and corresponding eigenvectors ${\mathbf{v}}_i={(v_{1,i},v_{2,i},\dots ,v_{n+1,i})}^T$ of the discrete Laplacian operator A (see Equation (10)) on $[0,L]$, coupled with homogeneous Neumann boundary conditions, obtained with a finite difference approximation on a uniform mesh of $n+1$ nodes with step size $h=L/n$, are ${\mu }_i={\displaystyle \frac{4}{h^2}}sin{\left({\displaystyle \frac{(i-1)\pi }{2n}}\right)}^2$ and $v_{j,i}=\left\{\begin{array}{cc}{\displaystyle \frac{1}{\sqrt{n}}},\hfill & \text{if}j=1\hfill \\ {\displaystyle \frac{1}{\sqrt{n}}}cos\left({\displaystyle \frac{\pi (i-1)(j-0.5)}{n}}\right),\hfill & \text{if}j=1,2,\dots ,n+1\hfill \end{array}\right.$ for $i=1,2,\dots ,n+1$.

Eigenvalues of the discrete fractional Laplacian operator, $A^{\gamma /2}$, can be taken as ${\left\{{\mu }_i^{\gamma }\right\}}_{i=1}^{n+1}$, since the Neumann boundary conditions of the discrete Laplacian, operator, A, properly transfer to the discrete fractional Laplacian operator, $A^{\gamma /2}$. Another representation for the fractional order Laplace operator on $[0,L]$ subject to standard Neumann Boundary conditions (the reflection boundary condition) has been presented in Cusimano et al. (2017) and this representation is called Reflection Matrix. The Reflection matrix of fractional order $\gamma$ on a uniform mesh of $n+1$ nodes with step size $h=L/n$ on $[0,L]$ is denoted as $A_{\gamma ,h}$. The Eigenvalues ${\nu }_{i,\gamma }$ of $A_{\gamma ,h}$ are derived in Cusimano et al. (2017) of the form:$${\nu }_{i,\gamma }={\left({\displaystyle \frac{2}{h}}sin\left({\displaystyle \frac{(i-1)\pi }{2n}}\right)\right)}^{\gamma }{\displaystyle \frac{1}{cos(\pi \gamma /2)}}cos\left({\displaystyle \frac{2(i-1)+\gamma (n+1-i)\pi }{2n}}\right),i=1,2,\dots ,n+1.$$

Then for $1<\gamma \le 2$, we can show that, as h goes to zero, both ${\left\{{\mu }_i^{\gamma }\right\}}_{i=1}^{n+1}$ and ${\left\{{\nu }_{i,\gamma }\right\}}_{i=1}^{n+1}$ converge to ${\left\{{\lambda }_i^{\gamma }\right\}}_{i=1}^{\infty }={\left\{{\left((i-1)\pi /L\right)}^{\gamma }\right\}}_{i=1}^{\infty }$, which are eigenvalues of the continuous Laplacian operator on $[0,L]$ subject to standard Neumann boundary conditions.

The discrete Laplacian on $[0,L_1]\times [0,L_2]$ can be represented in terms of discrete Laplacians on $[0,L_1]$ and $[0,L_2]$. Let $A_1$, $A_2$ and A be the discrete Laplacian operators on $[0,L_1]$, $[0,L_2]$ and $[0,L_1]\times [0,L_2]$, respectively. Boundary conditions embedded on $A_1$, $A_2$ and A should coincide. Then $A=A_1\otimes {\mathbf{I}}_y+{\mathbf{I}}_x\otimes A_2$, where ${\mathbf{I}}_x$ and ${\mathbf{I}}_y$ are the identity operators on $[0,L_1]$ and $[0,L_2]$, respectively. Based on this result we can obtain eigenvalues and eigenvectors of discrete Laplacian on $[0,L_1]\times [0,L_2]$. Let ${\lambda }_{\mathit{xy}}$ and ${\mathbf{v}}_{\mathit{xy}}$ denote the eigenvalues and eigenvectors of the discrete Laplace operator on $[0,L_1]\times [0,L_2]$ subject to homogeneous Neumann boundary conditions. Then we get ${\lambda }_{\mathit{xy}}={\lambda }_x\otimes {\mathbf{I}}_y+{\mathbf{I}}_x\otimes {\lambda }_y$ and ${\mathbf{v}}_{\mathit{xy}}={\mathbf{v}}_x\otimes {\mathbf{v}}_y$. Here ${\lambda }_x$, ${\lambda }_y$ are eigenvalues of the discrete Laplacian operators subject to homogeneous Neumann boundary conditions on $[0,L_1]$ and $[0,L_2]$ respectively. ${\mathbf{v}}_x$ and ${\mathbf{v}}_y$ are eigenvectors of the discrete Laplace operator subject to homogeneous Neumann boundary conditions on $[0,L_1]$ and $[0,L_2]$, respectively.

Consider the system of fractional diffusion Equations (11). Discretising first two equations of this system in time and treating fixed point iteration for each time step n we get the iteration processes:(21) $$\begin{array}{cc}\hfill {\displaystyle \frac{u^{n+1,m}-u^n}{\mathrm{\Delta }t}}& =-d_1{(-\mathrm{\Delta })}^{\gamma /2}{u}^{n+1,m}+F(u^{n+1,m-1},v^{n+1,m-1}),\hfill \\ \hfill {\displaystyle \frac{v^{n+1,m}-v^n}{\mathrm{\Delta }t}}& =-d_2{(-\mathrm{\Delta })}^{\gamma /2}{v}^{n+1,m}+G(u^{n+1,m-1},v^{n+1,m-1}).\hfill \end{array}$$(21)

Then the ith Fourier mode of the time-space discretisation of (21) becomes:(22) $$\begin{array}{cc}\hfill {\widehat{u}}^{n+1,m}& ={\displaystyle \frac{1}{1+d_1{\mu }_i^{\gamma /2}\mathrm{\Delta }t}}\left[{\widehat{u}}^n+\mathrm{\Delta }t\widehat{F}(u^{n+1,m-1},v^{n+1,m-1})\right],\hfill \\ \hfill {\widehat{v}}^{n+1,m}& ={\displaystyle \frac{1}{1+d_2{\mu }_i^{\gamma /2}\mathrm{\Delta }t}}\left[{\widehat{v}}^n+\mathrm{\Delta }t\widehat{G}(u^{n+1,m-1},v^{n+1,m-1})\right],\hfill \end{array}$$(22)

where $\widehat{u}$, $\widehat{v}$, $\widehat{F}$ and $\widehat{G}$ denote the Fourier coefficients of u, v, F and G, respectively, and ${\mu }_i$ is the ith eigenvalue of the discrete Laplace operator subject to considered boundary conditions.

At each time step, iterative process (22) is initiated taking ${\widehat{u}}^{n+1,0}={\widehat{u}}^n$ and ${\widehat{v}}^{n+1,0}={\widehat{v}}^n$. After calculating ${\widehat{u}}^n$ and ${\widehat{v}}^n$, the inverse reconstructions $u^n$ and $v^n$ in physical space can be computed by inverse Fourier transforms for the one-dimensional and two-dimensional cases.

4.2. Some other possible numerical methods

Discretising in space, a fractional reaction diffusion equation can be approximated to a system of ODEs using Matrix Transfer Technique. Then the system of ODEs can be solved using suitable numerical methods. In this section we briefly explain three possible numerical methods.

Applying the Matrix Transfer Technique, the first two equations of the system (11) can be presented in the form(23) $$\begin{array}{cc}\hfill {\displaystyle \frac{\text{d}\mathbf{u}}{\text{d}t}}& =-A^{\gamma /2}\mathbf{u}+\mathbf{F}(\mathbf{u},\mathbf{v}),\hfill \\ \hfill {\displaystyle \frac{\text{d}\mathbf{v}}{\text{d}t}}& =-d{A}^{\gamma /2}\mathbf{v}+\mathbf{G}(\mathbf{u},\mathbf{v}).\hfill \end{array}$$(23)

Here A is the discrete Laplacian operator coupled with standard boundary conditions. In one dimension, subject to Neumann boundary conditions, the matrix A is given in Equation (10). The vectors $\mathbf{u}$, $\mathbf{v}$, $\mathbf{F}$ and $\mathbf{G}$ represent the spatial discretisations of u, v, F and G, respectively. $A^{\gamma /2}=V{D}^{\gamma /2}{V}^{-1}$, where V is the matrix in which columns are eigenvectors of A, and D is a diagonal matrix in which elements are the eigenvalues of A.

4.2.1. A semi implicit scheme

Consider the following time discretisation of the system (23):(24) $$\begin{array}{cc}\hfill {\displaystyle \frac{{\mathbf{u}}^{n+1}-{\mathbf{u}}^n}{\mathrm{\Delta }t}}& =-A^{\gamma /2}{\mathbf{u}}^{n+1}+\mathbf{F}({\mathbf{u}}^n,{\mathbf{v}}^n),\hfill \\ \hfill {\displaystyle \frac{{\mathbf{v}}^{n+1}-{\mathbf{v}}^n}{\mathrm{\Delta }t}}& =-d{A}^{\gamma /2}{\mathbf{v}}^n+\mathbf{G}({\mathbf{u}}^n,{\mathbf{v}}^n).\hfill \end{array}$$(24)

This system can be arrange in the form:(25) $$\begin{array}{cc}\hfill (I+\mathrm{\Delta }t{A}^{\gamma /2}){\mathbf{u}}^{n+1}& ={\mathbf{u}}^n+\mathrm{\Delta }t\mathbf{F}({\mathbf{u}}^n,{\mathbf{v}}^n),\hfill \\ \hfill (I+\mathrm{\Delta }td{A}^{\gamma /2}){\mathbf{v}}^{n+1}& ={\mathbf{v}}^n+\mathrm{\Delta }t\mathbf{G}({\mathbf{u}}^n,{\mathbf{v}}^n),\hfill \end{array}$$(25)

where $A^{\gamma /2}=V{D}^{\gamma /2}{V}^{-1}$. This is a semi-implicit scheme and the systems of linear equations (25) were solved in MATLAB, using the Conjugate Gradient Scheme (CGS).

4.2.2. A fully explicit scheme

Substituting $A^{\gamma /2}=V{D}^{\gamma /2}{V}^{-1}$ in (25) we get:(26) $$\begin{array}{cc}\hfill (I+\mathrm{\Delta }tV{D}^{\gamma /2}{V}^{-1}){\mathbf{u}}^{n+1}& ={\mathbf{u}}^n+\mathrm{\Delta }t\mathbf{F}({\mathbf{u}}^n,{\mathbf{v}}^n),\hfill \\ \hfill (I+\mathrm{\Delta }tdV{D}^{\gamma /2}{V}^{-1}){\mathbf{v}}^{n+1}& ={\mathbf{v}}^n+\mathrm{\Delta }t\mathbf{G}({\mathbf{u}}^n,{\mathbf{v}}^n).\hfill \end{array}$$(26)

These systems of equations can be arranged in the form:(27) $$\begin{array}{cc}\hfill {\mathbf{w}}_1^{n+1}& =D_1^{-1}({\mathbf{w}}_1^n+\mathrm{\Delta }t{V}^{-1}\mathbf{F}(V{\mathbf{w}}_1^n,V{\mathbf{w}}_2^n)),\hfill \\ \hfill {\mathbf{w}}_2^{n+1}& =D_2^{-1}({\mathbf{w}}_2^n+\mathrm{\Delta }t{V}^{-1}\mathbf{G}(V{\mathbf{w}}_1^n,V{\mathbf{w}}_2^n)).\hfill \end{array}$$(27)

Where $D_1=(I+\mathrm{\Delta }t{D}^{\gamma /2})$, $D_2=(I+\mathrm{\Delta }td{D}^{\gamma /2})$ are diagonal matrices, ${\mathbf{w}}_1^n=V^{-1}{\mathbf{u}}^n$ and ${\mathbf{w}}_2^n=V^{-1}{\mathbf{v}}^n$. This is a fully explicit scheme. The system (27) can be solved for ${\mathbf{w}}_1^n$ and ${\mathbf{w}}_2^n$. ${\mathbf{u}}^n$ and ${\mathbf{v}}^n$ are given by ${\mathbf{u}}^n=V{\mathbf{w}}_1^n$, ${\mathbf{v}}^n=V{\mathbf{w}}_2^n$.

Also, the system of ODEs (23) was solved using the MATLAB function ODE15s that is suitable to solve stiff ODEs.

5. Numerical results

5.1. Fourier spectral method

We solved the system (8) numerically, on bounded one and two dimensional domains, subject to homogeneous Neumann boundary conditions with suitably chosen initial conditions using the numerical scheme (23). In these simulations, in each time step fixed point iterations proceed until the Relative Error (RE) becomes smaller than the predefined tolerance (Tol). The relative Error at each iteration m at time step n is calculated as follows:(28) $$\begin{array}{c}\hfill {\text{RE}}^{n,m}=\text{max}\left\{{\displaystyle \frac{\Vert {\mathbf{u}}^{n,m}-{\mathbf{u}}^{n,m-1}\Vert }{\Vert {\mathbf{u}}^{n-1}\Vert }},{\displaystyle \frac{\Vert {\mathbf{v}}^{n,m}-{\mathbf{v}}^{n,m-1}\Vert }{\Vert {\mathbf{v}}^{n-1}\Vert }}\right\}.\end{array}$$(28)

In the following simulations, the tolerance is set as Tol $={10}^{-13}$.

5.1.1. In one dimension

The system (22) was simulated on $[0,L]\subset \mathbb{R}$ in the time interval $[0,T]$ taking the following initial conditions(29) $$\begin{array}{c}\hfill u(x,0)=u_{02}x\in \mathrm{\Omega },\\ \hfill v(x,0)=\left\{\begin{array}{cc}{v}_{02}+v_{02}{\displaystyle \frac{\text{Rand}[0,1]}{10}},\hfill & x\in \omega \hfill \\ 0,\hfill & x\in \mathrm{\Omega }\backslash \omega .\hfill \end{array}\right.\end{array}$$(29)

Here $\omega =\left({\displaystyle \frac{L}{2}}-{\displaystyle \frac{2L}{100}},{\displaystyle \frac{L}{2}}+{\displaystyle \frac{2L}{100}}\right)$.

Figure 6. Density plot of v(x, t) for different domain sizes.

We controlled the termination of the iterative scheme (22) by applying the tolerance, Tol, to the relative error. Density plots of numerically evaluated v(x, t) for different values of $\gamma$ are shown in Figure 6 for $\alpha =5.2$, $\lambda =2.5$, $d=0.01$. In these simulations $L=8$, $T=100$, and the spatial and time grid sizes are $\mathrm{\Delta }x=0.05$ and $\mathrm{\Delta }t=0.002$, respectively.

Figure 7 illustrates the variation of the number of iterations in each time step. According to Figure 1, we can observe that, when $\gamma$ decreases, the required number of iterations of the iterative process (22) increases in order to maintain the fixed predefined error tolerance. As a result of that computational time increases. Figure 8 shows how the error defined in Equation (28), is maintained by the tolerance at each iteration. We can observe that as time passes the required number of iterations to maintain error tolerance decreases. Table 1 shows the total number of iterations, and total computational time of the above simulations.

Figure 7. Variation of the number of iterations against number of time steps when the step size is 0.002.

Figure 8. The error with respect to number of time steps when the step size is 0.002.

Table 1. Total number of iterations and computation time

$\gamma$	Total number of iterations	Computation time (in seconds)
2	411,513	144.53
1.8	457,521	158.66
1.6	558,109	190.89
1.4	682,707	233.26

According to Table 1, the total number of iterations required to maintain the tolerance, and computational time, increases as $\gamma$ decreases.

5.1.2. In two dimension

We solved the system (22) numerically on the 2D domain, $\mathrm{\Omega }\equiv [0,L]\times [0,L]$, subject to zero-flux boundary conditions and initial conditions:(30) $$\begin{array}{cc}& u(x,y,0)=u_{02}(x,y)\in \mathrm{\Omega },\hfill \\ \hfill & v(x,y,0)=\left\{\begin{array}{cc}{v}_{02}+v_{02}{\displaystyle \frac{\text{Rand}[0,1]}{10}},\hfill & (x,y)\in \omega \hfill \\ 0,\hfill & (x,y)\in \mathrm{\Omega }\backslash \omega .\hfill \end{array}\right.\hfill \end{array}$$(30)

Here $\omega =\left\{(x,y)\in \mathrm{\Omega },{\displaystyle \frac{L}{2}}-{\displaystyle \frac{2L}{100}}<x<{\displaystyle \frac{L}{2}}+{\displaystyle \frac{2L}{100}},{\displaystyle \frac{L}{2}}-{\displaystyle \frac{2L}{100}}<y<{\displaystyle \frac{L}{2}}+{\displaystyle \frac{2L}{100}}\right\}.$

Isosurfaces of numerically evaluated $v(x,y,t)$ for different levels of $\gamma$ are shown in Figures 9 for $\alpha =2.89$, $\lambda =6.1$, $d=0.01$. In these simulations, $L=8$, $T=60$, and the spatial and time grid sizes are $\mathrm{\Delta }x=0.05$ and $\mathrm{\Delta }t=0.002$, respectively. According to this figure, we can see that the heterogeneity of the space patterns increases as $\gamma$ decreases.

Figure 9. Isosurfaces of numerical solutions $v(x,y,t)$ for the parameter values with $\lambda =3.0$, $\alpha =6.1$, $d=0.01$.

5.2. Comparison of numerical schemes

In addition to the Fourier Spectral method, we solved the system of ODEs (23) using the MATLAB function ODE15s and numerical schemes (25) and (27) in one dimensional space subject to the same boundary conditions and initial data used in numerical simulations of Section 5.1.1. Density plots of the numerical solutions, obtained using these numerical methods, are presented in Figure 10. The parameters used in these simulations are $\alpha =4.201$, $\lambda =2.1$, $\gamma =1.6$ and $d=0.01$. We solved the linear systems arising at each time step of the numerical scheme (25), using the CGS method. The iteration processes of the CGS were controlled by predefined tolerance, TOL=${10}^{-13}$. We define the relative residual, RelRes, of a numerical solution of the linear system $Ax=b$ as RelRes $={\displaystyle \frac{\Vert b-Ax\Vert }{\Vert b\Vert }}$. Figure 11 shows the RelRes corresponding to the last iteration of the iterative process at each time step of the scheme (25). We can observe that RelRes of this scheme at each time step has been properly controlled by the TOL (that is the Rel-Res corresponding to last iteration of the each time step is less than the TOL). Therefore, numerical scheme (25) is convergent, and so the numerical solutions of the system (25), obtained by the above schemes, converge to a solution of the system (23).

Figure 10. Density plot of numerical solutions v(x, t) obtained by different Numerical Methods.

Figure 11. The relative residual plot of the last iteration of CGS used in semi implicit method.

We observed that all the density plots in Figure 10 are almost the same. Therefore all four numerical methods, the Fourier Spectral method, MATLAB ODE15s, Semi-implicit method and Fully-Explicit method give almost the same solution of system (11).

6. Discussion and conclusions

There are three spatially homogeneous steady states $S_1$, $S_2$ and $S_3$ of the model and Turing type instability occurs at $S_2$ only. We performed Turing instability analysis about $S_2$ and observed that instability conditions do not depend on the fractional power $\gamma$ but the dispersion relation and the growth rate do. Figure 2 shows the variation of Re($\sigma$) against the wave mode k for different levels of $\gamma$ and d. The range of unstable wave modes (the range of k) increases as $\gamma$ decreases when d is fixed. The wave mode correspond to the maximum growth rate, ${\sigma }_{\text{max}}$, increases as $\gamma$ decreases. Also for a fixed d, ${\sigma }_{\text{max}}$ does not vary with $\gamma$.

Figure 3 shows the variation of the real part of $\sigma$, Re($\sigma$), against k for different levels of d. According to this figure, we can observe that when $d>d_c$, Re($\sigma$) increases as d decreases from $d_c$. Also, the wave mode corresponds to the maximum growth rate and the range of unstable wave modes increases as d decreases from $d_c$. These observations imply that as d decreases from $d_c$, the growth rate, of coral increases, and the heterogeneity of the coral structures increases.

Figure 4 shows the variation of ${\sigma }_{\text{max}}$ with respect to $\alpha$ for $\lambda =2.1$. For fixed d and $\lambda$ the growth rate maximizes at $\alpha ={\alpha }_{\text{min}}$ (if $\lambda >2$) or at $\alpha =2\lambda$ (if $2d<\lambda \le 2$). For fixed d, as $\alpha$ increases from its minimum value (${\alpha }_{}\text{min}$ or $2\lambda$) to a maximum value ${\alpha }_{\text{max}}$ in the Turing space, the maximum growth rate ${\sigma }_{\text{max}}$ decreases from its maximum value to 0. Furthermore, we can observe that as d decreases, the maximum growth rate, ${\sigma }_{\text{max}}$, increases when $\alpha$ and $\lambda$ are fixed.

Figure 5 shows the plot of $n_{\text{max}}$ (wave number corresponding to ${\sigma }_{\text{max}}$) against $\alpha$ for different levels of $\gamma$ for $\lambda =2.1$, $d=0.01$. We observed that $n_{\text{max}}$ increases when $\gamma$ decreases while other parameters are fixed. Also, $n_{\text{max}}$ increases as $\alpha$ increases from ${\alpha }_{\text{min}}$ (or $2\lambda$) to ${\alpha }_{\text{max}}$ while other parameters are fixed. This implies that more heterogeneous spatial patterns are generated by the solutions at larger $\alpha$ and smaller $\gamma$ values.

The parameters d, $\alpha$ and $\lambda$ of the proposed model are defined as $d=d_v/d_u$, ${\alpha }^2=b_1{u}_s^2/b$ and $\lambda =b_2/b_1$. Turing type patterns are possible if these parameters lie in the Turing instability region. If we assume that the growth factors of coral are controlled except the nutrient availability, then it is reasonable to assume that $d_u$, $d_v$, b, $b_1$ and $b_2$ are fixed. Thus $\lambda$ and d are fixed and $\alpha$ depends on the nutrient supplying rate, a, only. By adjusting the nutrient supplying rate we can adjust $\alpha$ such that all parameters lie in the Turing instability region. Then spatial temporal patterns of coral can be expected when a sufficient perturbation from the steady state $S_2$ is done, in order to make the system unstable.

According to the results shown in Figure 5, as $\alpha$ increases the heterogeneity of the branching structures increases when d, $\lambda$ and $\gamma$ are fixed . That is if the nutrient supply rate increases then the heterogeneity of the branching structures increases.

Figures 6 and 9 shows the effect of the fractional order, $\gamma$, for the spatial patterns generated by the numerical solutions of the model. These numerical solutions show that the spatial heterogeneity of the solutions increases as $\gamma$ decreases. This behaviour matches with some of the theoretically derived results in Section 3.5. According to Figure 10 we can conclude that numerical solutions obtained by four numerical methods explained in Section 4 are almost same.

According to Henry and Wearne (2000) spatial sub-diffusion ($\gamma <2$) terms arise when modelling processes with spatial constraints such as porosity and fractal lattices, and super-diffusion ($\gamma >2$) terms arise when turbulent or chaotic behaviours appear in the media of the model. Based on this explanation we can give following physical interpretations on the effect of the parameter $\gamma$ of our model.

• If the porosity of the media increases (that means $\gamma$ decreases from 2 to 1) then the heterogeneity of the patterns increases.

• When the coral density of the tank increases the porosity increases. In other words as time passes the coral density increases. As a result of this, porosity increases and the heterogeneity of the coral patterns increases. Therefore it would be interesting to consider time dependent fractional orders which will be considered in future work.

• As the turbulence of the media increases (so that $\gamma$ increases from 2) then the heterogeneity of the coral patterns decreases.

Future work could consider time fractional components as well as calibration against experimental data. Refinements of the model could consider other factors of coral growth (such as light intensity to the reef, pH value of water, sedimentation, salinity, depth to the reef) and properties of the flow (flow velocity and turbulence).

Acknowledgements

The authors would like to thank Professor Fawang Liu, School of Mathematical Sciences, Queensland University of Technology, Brisbane, Australia for giving suggestions to improve the quality of the paper. The first author greatly appreciate the support given by the School of Mathematical Sciences, Queensland University of Technology, Brisbane, Australia through a visiting fellowship position at School of Mathematical Sciences to carry out this research.

Additional information

Funding

The authors received no direct funding for this research.

Notes on contributors

Lekam Watte Somathilake

Lekam Watte Somathilake is currently a Senior Lecturer in Mathematics at the Department of Mathematics, Faculty of Science, University of Ruhun, Matara, Sri Lanka. He completed PhD and MPhil degrees in Mathematics at University of Ruhuna, Sri Lanka in 2012 and 2004 respectively. His research interests include Ordinary and Partial Differential Equations, Reaction Diffusion Pattern Formations, Numerical Analysis and Mathematical Biology.

Kevin Burrage

Kevin Burrage is currently a Professor of Computational Maths at the School of Mathematical Sciences, Faculty of Science and Engineering, Queensland University of Technology, Brisbane, Australia. He is also a Visiting Professor to the Department of Computer Science, University of Oxford, UK. His research interests include Computational Mathematics, Computational Biology and Systems Biology. This research was carried out while L.W. Somathilake was serving as a Visiting Fellow (from July 2016 to June 2017) under the supervision of Professor Kevin Burrage at the Queensland University of Technology, Brisbane, Australia.