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

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 function ODE15s. 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. Subjects: Mathematical Modeling; Non-Linear Systems; Mathematical Biology


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.

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, 2003aMerks, , 2003bMerks, Hoekstra, Kaandorp, & Sloot, 2003a, 2003bMistr & Bercovici, 2003). A reaction diffusion type mathematical model for growth of corals in a tank has been proposed by 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.
The general form of the one variable fractional reaction diffusion equation is: Here −D (−Δ) ∕2 u, D and f (u, v) represent a fractional diffusion process, diffusion rate and the reaction process, respectively. For = 2, anomalous diffusion reduces to normal diffusion. For < 2 ( > 2), the diffusion process is slower (faster) than normal diffusion, where it is called sub-diffusive (super-diffusive). Subdiffusion 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. (1) 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.

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): 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: 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: This system can be written in the form: . That is b(u − u s ) is the overall nutrient supply rate to the system.
Time rate change of Nutrient concentration in reaction process = Anomalous diffusion of u + Constant supply rate of u − wasting rate of u − Reactive loss of u . (2) (3) 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: which gives 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: Finally we get a set of fractional reaction diffusion equations There are six parameters in this model. In order to reduce the number of parameters, we use nondimensionalisation techniques.

Nondimensionalisation
Now we nondimensionalise the equations system (6) by applying the coordinate transformations For notational convenience, eliminating bars and stars we obtain

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 ( = 2). That is, the boundary conditions can be expressed as: where denotes the outwards unit normal vector to the boundary Ω.
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. (2005Ilić et al. ( , 2006 plays an important role.
Let A be the discrete Laplacian coupled with standard homogeneous Neumann boundary conditions and let A = VDV −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: According to Ilić et al. (2005Ilić et al. ( , 2006, the discrete fractional operator A ∕2 (1 < ≤ 2) is defined as A ∕2 = VD ∕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 ∕2 .
Very recently, it has been proved that Neumann boundary conditions embedded in A transfers properly to A ∕2 in the one dimensional case (Cusimano, Burrage, Turner, & Kay, 2017).

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).

Steady states
Now we represent the first two equations of the system (8) of the form: There are three homogeneous steady The Jacobian matrix evaluated at S i is The trace and determinant of A are Tr( respectively.

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 and . 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).
Here, f (u, v), and g (u, v) are functions of u and v, representing the reaction processes of the relevant model, and Ω is the domain of the considered model. The Turing instability conditions of this model take the form: where and s ≡ (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 (12) The matrix form of this system is: Since S 2 is a stable steady state, C � 1 ≡ a 11 a 22 − a 12 a 21 > 0, C � 2 ≡ a 11 + a 22 < 0 and hence g( ) < 0. Therefore, the real parts of the solutions for in equation (13)  In order for spatial patterns to exist, h min < 0 and k c > 0. That is the conditions: 2d .
Therefore, the instability conditions for the fractional reaction diffusion equations, and for the standard Reaction Diffusion equations are same. But the growth rate, , and wave numbers, k, depend on . In terms of the model parameters, we can represent the instability conditions at S 2 in the following form:

Instability regions
It is clear that C 1 holds at S 2 in the region > 2 . The bifurcation C 2 = 0 (Hopf Bifurcation), gives one positive real roots for when > 1. These are Equation C 3 = 0 gives two real roots for when > d. These are Equation C 4 = 0 gives two positive real roots for when > d. These are Since > 2 , letting = 2 + ( > 0) in C 2 , we obtain Therefore C 2 is negative for ≤ 2. Also, (3) 1 > (4) 1 for ≠ 2d and (3) 1 = (4) 1 for = 2d. Moreover Therefore, C 3 is not satisfied in the region ≤ 2d. Also, we can show numerically that . (4) By considering these facts, we can write the instability region (Turing space) as a union of two regions  1 and  2 . Here, where min ( ) = (2) 1 and max ( , d) = (4) 1 . Figure 1 shows the parameter regions  1 and  2 in , , d space as volumes.

Critical values of d
√ 2) and > 2, the instability parameter region  2 is unbounded and when d > d c 1 and > 2 it is bounded (Somathilake & Wedagedera, 2014). In order to satisfy Turing instability conditions the diffusion rate d should satisfy d < d c ( , ) (Somathilake & Wedagedera, 2014).

Effect of the fractional exponent and model parameters for the growth
where We can show that is maximized at k = k 1 . That is max ( , ) = | k=k 1 . The plots of max against when = 2.1 for different levels of d are shown in Figure 4.
It is clear that max ( , ) is maximized at = min ( ) for given . The wavelength corresponding to the maximum growth rate is given by max = 2 ∕k 1 . The wave number n max corresponding to the maximum growth rate is given by n max = k 1 L∕ . The plot of n max against for = 2.1, d = 0.01, at  different levels of is shown in Figure 5 . We can observe that for fixed and d, the value of n max increases when increases from 1. That is the heterogeneity of the spatial patterns increases as increases when , and d are fixed. Also we observe that when both and increase the possible minimum and maximum values of n max as well as the range of n max increase. This implies that the heterogeneity of the spatial patterns increases as and increase when d and are fixed.

Fourier spectral method
The Fourier spectral method has been used in numerical computations of fractional reaction diffusion equations by several researchers (  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. (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

The eigenvalues i and corresponding eigenvectors
Eigenvalues of the discrete fractional Laplacian operator, A ∕2 , can be taken as i n+1 i=1 , since the Neumann boundary conditions of the discrete Laplacian, operator, A, properly transfer to the discrete fractional Laplacian operator, A ∕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 on a uniform mesh of n + 1 nodes with step size h = L∕n on [0, L] is denoted as A , h . The Eigenvalues i, of A , h are derived in Cusimano et al. (2017) of the form: Then for 1 < ≤ 2, we can show that, as h goes to zero, both i  Note: The values of other parameters are = 2.1, d = 0.01, 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: Then the ith Fourier mode of the time-space discretisation of (21) becomes: where û, v, F and Ĝ denote the Fourier coefficients of u, v, F and G, respectively, and 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 û n+1,0 =û n and v n+1,0 =v n . After calculating û n and 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.

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 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 , , and represent the spatial discretisations of u, v, F and G, respectively 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.

A semi implicit scheme
Consider the following time discretisation of the system (23): This system can be arrange in the form: where A ∕2 = VD ∕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).

A fully explicit scheme
Substituting A ∕2 = VD ∕2 V −1 in (25) we get: These systems of equations can be arranged in the form: Where D 1 = (I + ΔtD ∕2 ), D 2 = (I + ΔtdD ∕2 ) are diagonal matrices, n 1 = V −1 n and n 2 = V −1 n . This is a fully explicit scheme. The system (27) can be solved for n 1 and n 2 . n and n are given by Also, the system of ODEs (23) was solved using the MATLAB function ODE15s that is suitable to solve stiff ODEs.

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: In the following simulations, the tolerance is set as Tol = 10 −13 .

In one dimension
The system (22) Figure 7, we can observe that, when 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 (26) (I + ΔtVD ∕2 V −1 ) n+1 = n + Δt ( n , n ), maintain error tolerance decreases. Table 1 shows the total number of iterations, and total computational time of the above simulations.
According to Table 1, the total number of iterations required to maintain the tolerance, and computational time, increases as decreases.   Isosurfaces of numerically evaluated v(x, y, t) for different levels of are shown in Figures 9 for = 2.89, = 6.1, d = 0.01. In these simulations, L = 8, T = 60, and the spatial and time grid sizes are Δx = 0.05 and Δt = 0.002, respectively. According to this figure, we can see that the heterogeneity of the space patterns increases as decreases.

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 = 4.201, = 2.1, = 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 = ‖b − Ax‖ ‖b‖ . 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). 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).

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 but the dispersion relation and the growth rate do. Figure 2 shows the variation of Re( ) against the wave mode k for different levels of and d. The range of unstable wave modes (the range of k) increases as decreases when d is fixed. The wave mode correspond to the maximum growth rate, max , increases as decreases. Also for a fixed d, max does not vary with . Figure 3 shows the variation of the real part of , Re( ), against k for different levels of d. According to this figure, we can observe that when d > d c , Re( ) 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 max with respect to for = 2.1. For fixed d and the growth rate maximizes at = min (if > 2) or at = 2 (if 2d < ≤ 2). For fixed d, as increases from its minimum value ( min or 2 ) to a maximum value max in the Turing space, the maximum growth rate max decreases from its maximum value to 0. Furthermore, we can observe that as d decreases, the maximum growth rate, max , increases when and are fixed. Figure 5 shows the plot of n max (wave number corresponding to max ) against for different levels of for = 2.1, d = 0.01. We observed that n max increases when decreases while other parameters are fixed. Also, n max increases as increases from min (or 2 ) to max while other parameters are fixed. This implies that more heterogeneous spatial patterns are generated by the solutions at larger and smaller values.
The parameters d, and of the proposed model are defined as d = d v ∕d u , 2 = b 1 u 2 s ∕b and = 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 and d are fixed and depends on the nutrient supplying rate, a, only. By adjusting the nutrient supplying rate we can adjust 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 increases the heterogeneity of the branching structures increases when d, and 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, , 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 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 ( < 2) terms arise when modelling processes with spatial constraints such as porosity and fractal lattices, and super-diffusion ( > 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 of our model.
• If the porosity of the media increases (that means 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 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).