Dynamics of the time-fractional reaction–diffusion coupled equations in biological and chemical processes

This paper aims to demonstrate a numerical strategy via finite difference formulations for time fractional reaction–diffusion models which are ubiquitous in chemical and biological phenomena. The time-fractional derivative is considered in the Caputo sense for both linear and nonlinear problems. First, the Caputo derivative is replaced with a quadrature formula, then an implicit method is used for the remaining part. In the linear case, the proposed strategy reduces the time fractional models into linear simultaneous equations. In nonlinear cases, Quasilinearization is utilized to tackle the nonlinear parts. With this strategy, solutions of the fractional system transform into linear algebraic systems which are easy to solve. Next, the Von Neumann method is implemented to examine the stability of the scheme which discloses that the scheme is unconditionally stable. Further, the applicability of the presented scheme is tested with different linear and nonlinear models which include the one dimensional Schnakenberg and Gray–Scott models, and one and two dimensional Brusselator models. To analyze the accuracy of the present technique two norms namely, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathbb {L}_{\infty }$$\end{document}L∞ and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathbb {L}_{2}$$\end{document}L2, and relative error are addressed. Moreover, the obtained outcomes are shown tabulated and graphically which identifies that the scheme properly works for the time fractional reaction–diffusion systems.


Basic definitions
In this section of the manuscript, some basic definitions related to FC are discussed.

Definition
The fractional derivative of χ in Caputo sense 37 is defined as follow:

Definition
The one and two parameters Mittag-Leffler function is defined as follows 37 : The Caputo derivative of exponential and trigonometric functions can be expressed in the form of the Mittag-Leffler function which is described below: , where n = ⌈α⌉ , where n = ⌈α⌉.

Description of the method for one-dimensional TFRDMs
Consider the following TFRDM [38][39][40][41] : In Eq. ( 5) U (x, t) and V (x, t) are two dependent variables describes the dynamics where 0 < α ≤ 1 , a j , b j , c j , d j , e j , m j and κ j , are real constants for each j = 1, 2 which described some physical interpretation discussed in problem 5.1.F 1 and F 2 are the source terms and = [x 0 , x N ] is the spacial domain is divided into M subintervals where width of each subinterval is dx = x N −x 0 M .The associated initial conditions (ICs) are: The corresponding boundary conditions (BCs) are categorized in the following way: • Type 1: • Type 2: • Type 3: Now, based on the variation of the coefficient, Eq. ( 5) reduces to a variety of different linear and non-linear models.The models which will be under investigation, are listed in Table 1 [38][39][40][41] : Using Eq. ( 2) and an implicit scheme to Eq. ( 5) the resultant is: ( (2) (5) where Further simplification of the above Eqs.(10-11) leads to: The nonlinear terms , and U V 2 n+1 are linearized using Qausilinearization 42 : Inserting the central difference approximation of U n+1 xx and V n+1 xx together with non-linear terms Eqs.(14-16), in Eqs. ( 12) and ( 13) eventually gives linear system of equations with the compact form given by: where the dimensions of the matrices of A and W are described below: • For type 1 boundary conditions, the order of A and W will be respectively.• For type 2 boundary conditions, the order of A and W will be 2M × 2M and 2M × 1 , respectively.

Stability analysis
In this part of the manuscript, the stability analysis of the numerical method for the TFRDMs model is discussed.Assume the coefficients for this model from the first row of Table 1, then the system reduces to: (10)

Proof
Following the procedure 43 we assume: where i = √ −1. and some algebraic manipulation leads to: From Eq. ( 22) we get: For n = 0 In a similar fashion the following hold: If two consecutive approximations are closed then their difference approaches zero hence: Continuing in this way one can write Now putting Eq. ( 23) in Eq. ( 21) we obtain: which gives: (18)

Algorithm
Input: 0 < α ≤ 1 , a j , b j , c j , d j , e j , m j and κ j for j=1,2 and time step size τ.
Output:Solve the system of PDEs numerically using FDM to evaluate approximate solutions U n and V n .
Step 1: Generate computational grid over the domain , discretizing the spatial dimensions into a finite set of points.
Step 2: Apply central differences formula to discretize the spatial derivatives in the PDEs and quadrature formula for time derivative.
Step 4: Calculate matrices A and W n .
Step 5:Calculate W n+1 by using inversion method from A W n+1 = W n ,.

Numerical experiments
In this section, the proposed scheme is applied to some linear and non-linear RDMs.To assure the accuracy, global relative error (RE), L 2 and L ∞ have been used which are defined below: where χ n+1 j and χ n j are the approximate solutions at two consecutive time levels.
where χ ext j and χ app j are the exact and approximate solutions, respectively.

Problem 5.1 (linear model)
Consider the coefficients from the first row of Table 1 the resultant equation is given by [38][39][40][41] : The corresponding ICs and BCs are: The closed-form solutions of this system are: Using the fact ( ) the associated source terms are extracted as follows: (28) www.nature.com/scientificreports/For numerical simulations three cases are addressed here.

Diffusion-dominated
For diffusion-dominated case the parameters considered are a = 0.1, b = 0.01 and d = 1.

Reaction-dominated
For the reaction-dominated case the parameters considered are a = 2, b = 1 and d = 0.001

Reaction-dominated with stiff reaction
For this case the selection of parameters are a = 100, b = 1 and d = 0.001.This problem has been solved numerically, with the suggested technique, and the obtained results are noted in the form of tabulated and graphical forms.In Tables 2, 3, 4, 5, 6 and 7, the L 2 and L ∞ error norms are recorded for different times and α .The consecutive Tables 2 and 3 show the results for diffusion-dominated case, Tables 4  and 5 for reaction-dominated case and Tables 6 and 7 for reaction-dominated with stiff reaction case, respectively.Tabulated simulations reveal the good performance of the present technique.Similarly, solutions profile of exact versus numerical solutions are displayed in Figs. 1, 2 and 3 for diffusion dominated, reaction dominated and

Problem 5.2 (Brusselator model)
The exact solution this model and the next three nonlinear models is not available.For these cases ignored the corresponding source term.The numerical results in non-fractional form are given in literature.We tried the proposed technique for different values of α and observed its graphical behavior which approaches towards the integer value of α = 1 .In these cases the relative error is measured between two consecutive time levels.On behalf of these arguments we claim that proposed method works for such nonlinear fractional models.Using the coefficients from the second row of Table 1, which gives the following Brusselator model of kinetics for two chemical components [38][39][40][41] : with ICs: and the BCs are: (36) Table 5. L 2 norms for distinct values of α at different time levels in the reaction-dominated when when a = 2, b = 1, d = 0.001.8. Similarly, the density values for periodical motion are given in Table 9. Tabulated data discloses that the proposed technique gives good results in terms of error values which are comparable with the previous  Here, we take the values of coefficients from the third row of Table 1 which gives the Gray-Scot model [38][39][40][41] 10.Tabulated simulation reveals that the scheme works well for large time.Numerical solutions are plotted in Fig. 6 for integer cases which show that the outcome is in good agreement with available results in previous work.Likewise, in Fig. 7 the simulations are noted for fractional cases which indicate that the graphical behavior is nearly similar to integer when α = 0.9 .In Fig. 8 three dimensions numerical solutions are shown which show the wave type pattern.

Problem 5.4 (Schnakenberg model)
Choosing the coefficient values from the fourth row of Table 1, we obtain the following Schnakenberg model [38][39][40][41] : where U and V denote the activator and inhibitor concentrations, respectively, and d is diffusion coefficients, γ , a, and b are biological reaction rate constants.The following are the associated ICs:    Using the stratagem discussed earlier, Eq. ( 45) reduces to: (47) where Plugging the quadrature rule for fractional derivative in Eq. ( 48) and Eq. ( 49) the new equations are given below: Numerical illustration for Schnakenberg model in one dimensional graph for U and V for γ = 5000 and γ = 10,000 at t=2.5 when dt = 0.001 , dx = 0.005 and α = 0.75, 0.9 and 1.
Vol Rearranging the terms in the above system the results are: Linearization of the nonlinear term U 2 V n+1 is tackled by the following formula 42 : Inserting the approximation of the nonlinear term, and the central difference approximation for the involved derivative, the linear system of equations can be obtained which are given in compact form: where (51) (52) (54) Here, the following two-dimensional Brusselator model of kinetics for two chemical components is considered for validation of the scheme 44,45 : along with the following conditions: The associated source terms can be adjusted using the exact solution and the formula: where E 1,n−α+1 (:) is the Mittag-Leffler function defined earlier.The closed-form solution to the above problem is: We solve this model for the parameters ε 1 = ε 2 = 0.25, A = 0, B = 1, dx = dt .Also, the spatial and temporal domains are [0, 1] × [0, 1], and [0,5], respectively.In Table 12 the obtained error norms are matched with the previous norms given in the papers for integer case 44,45 .It is noticed that computed outcomes are matchable with available solutions.the scheme is tested for fractional values of α and the achieved results are reported in the form of L ∞ and L 2 norms in Tables 13 and 14, respectively.From these tables one can see, that the proposed scheme works for fractional cases as well.For further clarification, the solutions are sketched for integer and fractional cases in Figs.11 and 12, which show that exact and numerical solutions are promised well.

Conclusions and future plan
In this work, an implicit scheme has been addressed to solve RDCMs in fractional form.The involved fraction derivative and spatial derivatives were approximated with a well-known L 1 formula (Quadrature rule) and finite differences, respectively.Next, the stability of the scheme was investigated via Von Neumann analysis.Moreover, the scheme has been tested with different linear and nonlinear problems and the outcomes were compared with exact and existing results in literature.From tabulated simulations and graphical solutions, it has been observed that the proposed scheme works well for RDMs and can be used for such complicated problems having no exact solutions.In future the proposed methodology can be extended to the three dimensional problems coupling with different variable order fractional derivatives like, Caputo Fabrizio, and Atangana-Baleanu-Caputo etc.Moreover, the strategy can be tested for variable order local fractional derivative problems as well. (57) ∂ α U (x, y, t) ∂t α = ε 1 U xx + U yy − B + 1 U + U 2 V + A + F 1 (x, y, t), ∂ α V (x, y, t) ∂t α = ε 2 V xx + V yy + B U − U 2 V + F 2 (x, y, t), (58) U (x, y, 0) = e (−x−y) V (x, y, 0) = e (x+y) , x, y ∈ �.

Figure 1 .
Figure 1.Solution profile with absolute error of U and V at t = 0.001 for diffusion dominant case.

Figure 2 .
Figure 2. Solutions profile of U and V at t = 0.001 for reaction dominant case.

Figure 3 .
Figure 3. Solutions profile of U and V at t = 0.001 for reaction dominant with stiff reaction.

Figure 4 .
Figure 4. Numerical illustration of Brusselator model for U and V when α = 1..

Figure 12 .
Figure 12.Solution profiles of U and V at t=0.1 when α = 0.5 and dt = 0.001.

Table 3 .
L 2 norms for distinct values of α at different time levels in the diffusion-dominated when a = 0.1, b = 0.01, d = 1.

Table 4
. L ∞ norms for distinct values of α at different time levels in the reaction-dominated when when a = 2, b = 1, d = 0.001.

Table 6 .
L ∞ norms for distinct values of α at different time levels in the reaction-dominated with stiff when when a = 100, b = 1, d = 0.001.

Table 8 .
Relative error for various α values at various time levels when dt = 0.01 and dx = 0.05 of the Brusselator model.

Table 12 .
L ∞ norms for U and V at t = 2 and dt = 0.001 of Brusselator model in 2D for α = 1..

Table 13 .
L ∞ norm of Brusselator model in 2D for different α and time.