Fractional Reaction-Diffusion Equations for Modelling Complex Biological Patterns

We consider a system of nonlinear time-fractional reaction-diffusion equations (TFRDE) on a finite spatial domain x ∈ [0, L], and time t ∈ [0, T]. The system of standard reaction-diffusion equations with Neumann boundary conditions are generalized by replacing the first-order time derivatives with Caputo time-fractional derivatives of order α ∈ (0, 1). We solve the TFRDE numerically using Grünwald-Letnikov derivative approximation for timefractional derivative and centred difference approximation for spatial derivative. We discuss the numerical results and propose the applications of TFRDE for the modelling of complex patterns in biological systems. | Fractional reaction-diffusion | Fractional derivative | Caputo derivative | Grünwald-Letnikov derivative | ® 2012 Ibnu Sina Institute. All rights reserved. http://dx.doi.org/10.11113/mjfas.v8n3.135


INTRODUCTION
In recent years, the system of fractional reactiondiffusion (FRD) equations have gained considerable popularity to study nonlinear phenomena arise in the disciplines of science and engineering [1].Of these particular interests are patterns formations [2][3][4][5][6].The evolution of pattern formation is best described by the fractional-order models because the fractional derivatives take into consideration the whole history of the system which is called the memory effect [7].
The fractional reaction-diffusion equation is a generalization of the standard reaction-diffusion equation with derivative of arbitrary real order.The fractional reaction-diffusion equation is obtained by replacing the first-order time derivative index by α ∈ (0,1), or the second-order spatial derivative index by β ∈ (1,2), or both in the standard reaction-diffusion equation, (1) where and are the fractional derivative operators, κ is the dimensionless diffusion coefficient and denotes reaction kinetics.

Corresponding author at: E-mail addresses: kuazlina@siswa.um.edu.my (Ku Azlina Ku Akil)
There are several definitions of the fractional derivatives in literature such as the Riemann-Liouville, the Grünwald-Letnikov and the Caputo derivatives [7,8].
In this paper we consider a o ne-dimensional system of nonlinear time-fractional reaction-diffusion equations (TFRDE) on a finite spatial domain and time with Neumann boundary conditions.We solve the system of TFRDE numerically using Grünwald-Letnikov derivative approximation for time-fractional derivative and centred difference approximation for spatial derivative.This paper is organized as follows.In Section 2, we present the mathematical model of this study.The numerical scheme will be outlined in Section 3. The computational results are discussed in Section 4 and the conclusion will be drawn in Section 5.

MATHEMATICAL MODEL
In this section we introduce a one-dimensional system of nonlinear time-fractional reaction-diffusion equations (TFRDE) in the following form, , . (3) The above system subjects to the zero-flux boundary conditions at both ends of the spatial domain, Here, and denote the concentration for the two species, and are time characteristics of the system, and and are in general the nonlinear functions defined by , ( 5) where β and η are external parameters.The fractional derivatives and are the Caputo fractional derivatives in time of order , which are defined as (7) for is a Caputo fractional derivative operator.It should be noted that when α = 1, the system of equations ( 2) and (3) correspond to the system of standard reaction-diffusion equations.We have chosen the fractional derivative in the Caputo sense because it allows the utilisation of the initial condition into the formulation of the problem [7], which is similar to the integer-order differential equations.

NUMERICAL SCHEME
This section describes the numerical techniques to solve the system of TFRDE (2) and (3).We employed the numerical scheme as described by Gafiychuk et al. [4] and Ciesielski and Leszczynski [16].The system of TFRDE ( 2) and (3) can be written as a single partial differential equation as (8) for The grid points in the space domain [0,L] and the time domain [0,T] are labeled and , respectively, where is the grid size in space and is the grid size in time.The notation represents the value of function u at the grid point .The second order spatial derivative was approximated using the centred difference scheme, .
The time-fractional derivative was discretized using the definition of Grünwald-Letnikov (GL) approximation due to the fractional derivative defined in the GL sense is more flexible and straightforward for the numerical calculation purposes.The Caputo time fractional derivative is expressed in terms of the Riemann-Liouville (RL) fractional derivative, , for where The RL fractional derivative is equivalent to the definition of GL fractional derivative.The GL fractional derivative is defined as (11) The fractional derivative (11) can be approximated within the interval with subinterval as The term is known as GL coefficient with the following properties and ( 13) Applying ( 9) -( 13), the numerical scheme for equation ( 8) becomes (14) where , for We developed a semi implicit scheme in which both the spatial and the time-fractional derivatives were discretised implicitly.The stability and the accuracy of the method have been discussed in [15].The numerical scheme for each time layer for equation ( 14) represents a system of algebraic equations with block diagonal matrix and was solved at each iterations using the Newton-Raphson method [24].

RESULTS & DISCUSSION
The computer code for numerical scheme obtained in Section 3 was written in MATLAB 7.9.0(R2009b).The time and the spatial steps used in the simulation were varying from 0.005 to 0.1 and from 0.01 to 0.1, respectively.The initial conditions are similar to [5], , .and are the steady-state solution of equations ( 2) and (3) which correspond to homogeneous equilibrium state .Figures 1 and 2 show the results of computer solution of TFRDE considered above for different values of fractional index, α∈ (0,1) on the spatial domain and time domain Here, the values of time characteristics are fixed at unity, τ 1 = τ 2 = 1.As the value of α increases from 0.1 to 0.9 with step 0.1, we can see the evolution of pattern formation in the steady-state solutions of and and the system is stable.Fig. 3 The dynamics of variables u 1 (x,t) (top panel) and u 2 (x,t) (bottom panel) when the ratio is 0.135 for α = 0.6 to α = 0.9.The parameters are β = 2.0, η = -0.1,κ 1 = 0.05, κ 2 = 1.0,Our results in Figures 1 and 2 s how that α ≤ 0.5 the amplitude of both variables u 1 and u 2 increases steadily.However when α > 0.5 the amplitude increases drastically within time interval [0,10] and later becomes monotonous as time evolves.Furthermore, this reveals that the history of the system which is called memory effect plays significant role in the formation of pattern.
Figure 3 illustrates the transition of pattern formation in variables u 1 and u 2 from the steady-state structures (figures (a) and (b)) to the homogeneous oscillatory structures (figures (c) and (d)) that eventually destroy the pattern formation when the value α increases from 0.6 to 0.9 and the ratio of time characteristics, was fixed at 0.135.This shows that for α < 1, when the ratio of time characteristics less than unity and sufficiently small, the system becomes unstable.These temporal patterns behaviour have deep physical meaning that need to be investigated further such as the amount of energy transfer due to nonlinear interaction at each layer of the system.

CONCLUSION
A system of TFRDE with cubic nonlinearity has been numerically studied on finite spatial and time domains.We have demonstrated the temporal patterns behaviour in the steady-state solutions of system of TFRDE when the fractional derivative index α increases from 0.1 to 0.9.We have also showed that at certain value of α the formation of patterns transform from the steady-state structure to the homogenous oscillatory structure, given that the ratio of time characteristics is sufficiently small.The system of fractional reaction-diffusion equations has proven useful in understanding the dynamics of nonlinear phenomena.