Turing Instability of Gray-Scott Reaction-Diffusion Model with Time Delay Effects

Suqi Ma

Department of Mathematics, China Agricultural University, Beijing, China.

DOI: 10.4236/ijmnta.2023.122004                

Abstract

The reaction diffusion Gray-Scott model with time delay is put forward with the assumption of Neumann boundary condition is satisfied. Based on the Turing bifurcation condition, the Turing curves on two parameter plane are discussed without time delay. The normal form is computed via applying Lyapunov-Schmidt reduction method in system of PDE, and the bifurcating direction of pitchfork bifurcation underlying codimension-1 singularity of Turing point is computed. The continuation of Pitchfork bifurcation is simulated with varying free parameter continuously near the turing point, which is in coincidence with the theoritical analysis results. The wave pattern formation in the case of turing instability is also simulated which discover Turing oscillation phenomena from periodicity to irregularity.

Keywords

Reaction Diffusion, Turing Bifurcation, Normal Form, Time Delay

1. Introduction

Reaction diffusion model is ubiquitous in describing the spatial-temporal dynamical evolutionary behavior and to discover a variety of wave patterns which as often are seen in biological species in real life. The Gray-Scott diffusion model is motivated to do such simulations and have attracted attention in many researchers investigation work. Initially it was set forth by Gray and Scott as a variant of the autocatalytic model of glycolysis proposed by Sel’kov. Later authors in the papers ( [1] [2] [3] [4] ) put forth its new reaction mechanism by an auto-catalytic sequence and modelling its state variable control by time delay feedback method. For simplicity, the Gray and Scott model comes from mathematically model governed by

$\begin{array}{l}\frac{\text{d}{u}_1}{\text{d}t}=\epsilon d\Delta u_1+a-u_1-\frac{4u_1\left(t-\tau \right)u_2\left(t-\tau \right)}{1+u_1{\left(t-\tau \right)}^2}\\ \frac{\text{d}{u}_2}{\text{d}t}=d\Delta u_2+b{u}_1-b\frac{u_1\left(t-\tau \right)u_2\left(t-\tau \right)}{1+u_1{\left(t-\tau \right)}^2}\end{array}$ (1.1)

with the initial conditions and Neumann boundary condition

$\begin{array}{l}\frac{\partial u_1}{\partial n}=\frac{\partial u_2}{\partial n}=0,x\in \partial \Omega ,t>0\\ u_1\left(x,\theta \right)={\varphi }_1,u_2\left(x,\theta \right)={\varphi }_2,x\in \Omega ,\theta \in \left[-\tau ,0\right]\end{array}$ (1.2)

where $a\mathrm{,}b\mathrm{,}d\mathrm{,}\epsilon$ are constants, $\Delta$ is Laplacian operator and $\Omega$ is the domain of space variable x.

It is easily to calculate that Equation (1.1) has a homogeneous stationary solution

$u_1=\frac{ab}{b+4\sigma },u_2=\frac{\sigma \left(a^2{b}^2+b^2+8b\sigma +16{\sigma }^2\right)}{b{\left(b+4\sigma \right)}^2}$ (1.3)

Motivated by the aims to discover the complex dynamical evolution behavior both of the inhomogeneous solutions, we investigate the turing bifurcation mechanism inherently with time delay effects. Lyapunov-Schmidt method can be applied to investigate the bifurcation behavior at the Turing-instability bifurcation point with simple eigenvalue. Some authors as referring to the papers ( [5] [6] [7] [8] ) also develop the numerical algorithm of Lyapunov-Schmidt method to explore the bifurcation scenario near simple bifurcation point. According to Lyapunov-Schmidt Reduction method, the original partial delay differential equation is expressed on the invariant center manifold to acquire the norm form correspondingly. With the analyzing results of the characteristic equation, for example the roots with zero real part crossing imaginary axis underlying the related positive or negative transversal conditions, the turing bifurcation mechanism is exploited ( [9] [10] [11] [12] ). The numerical simulation results verify the analyzing result and near the threshold value the simulation bifurcating solution is in coincidence with its bifurcating directions.

We get the periodical solutions in space near the turing bifurcation point and time-periodically solution mainly dependent on a series of DDEs by discretion method underlying time delay. The dynamics of Turing patterns in one dimension are simulated, which reflects the spatio-temporal oscillation under time delay feedback. The wave simulations of reaction diffusion equation with time delay still take some interesting methods, which alike the well known differential quadrature algorithms ( [5] [13] [14] [15] [16] ), element free Galerkin method, and Trigonometric B-spline functions interpolation method, etc. With free parameters varying in Turing instability region, the periodical travelling wave pattern is induced ( [17] [18] ), and even from periodicity to un-regularity of wave patterns are discovered.

The whole paper is arranged as the listed. In section 2, the turing instability of Gray and Scott model underlying time delay is discussed. In section 3, the pitchfork bifurcation branch of turing bifurcation is discussed, which also determine the bifurcation direction in turing point. In section 4, the numerical simulation is done, the results is in coincidence with the theoretical analysis proof. Finally the discussion is given briefly.

2. Turing Instability

To solve Turing bifurcation problem, we first explains some notations in diffusion Equation (1.2). Suppose X is the domain of Laplacian operator and respectively, $X\mathrm{,}Y$ is Hilbert space defined by the inner product $\left(\varphi \mathrm{,}\psi \right)={\displaystyle {\int }_{\Omega }}\varphi {\left(\pi x\right)}^{\text{T}}\psi \left(\pi x\right)\text{d}x$ wherein $\varphi$ is defined in $\Omega ={\left[\mathrm{0,1}\right]}^n$ with $n\le 3$ . Mapping $G_{\mathrm{1,2}}\mathrm{<mo>\; :\; </mo>}\left(X\mathrm{,}\alpha \right)\to Y$ with $\alpha \in R^2$ , we write the elliptic Equation of (1.2) as

$\begin{array}{l}{G}_1\left(u,v\right)=\epsilon d\Delta u_1+a-u_1-\frac{4u_1{u}_2}{1+u_1^2}\\ G_2\left(u,v\right)=d\Delta u_2+b{u}_1-b\frac{u_1{u}_2}{1+u_1^2}\end{array}$ (2.1)

with $\left(u\mathrm{,}v\right)\in X$ and $\alpha =\left(\epsilon \mathrm{,}d\right)$ .

Notice we set time delay $\tau =0$ in Equation (2.1), the discussion of Turing instability of Equation (2.1) is independent of time delay. The turing-Hopf instability induce the periodical spatial oscillation phenomena doesn’t discuss here, however it may happen with time delay varying due to complex wave pattern phenomena.

To analyze turing bifurcation point, we firstly investigate the characteristic equation with simple zero root. Based on the known knowledge of diffusion equation, the characteristic equation of Equation (2.1) is written as ${\displaystyle {\prod }_{k=0}^N{W}_k\left(\lambda \right)}=0$ , with different $k\in N_0$ . It is easily to compute that

$\begin{array}{c}{W}_k\left(\lambda \right)=\left(\left(-4d{k}^2{\pi }^2-4\lambda \right)a^2+5b\left(d\epsilon k^2{\pi }^2+\lambda +5\right)a+100d{k}^2{\pi }^2+100\lambda \right)\\ +\left(a^2+25\right)\left(d{k}^2{\pi }^2+\lambda \right)\left(d\epsilon k^2{\pi }^2+\lambda +1\right)\\ =0\end{array}$ (2.2)

Turing instability expands the non-homogeneous solution bifurcation branches with spectral assumption at singularity point. The characteristic equation with zero root with $k=0$ is called as long wave modulation, however the named Turing instability happens if zero root appears with $k\ge 1$ . We set

$\begin{array}{l}Tr=\frac{a^{2}d\epsilon k^2+a^{2}d{k}^2+25d\epsilon k^2+5ab+25d{k}^2-3a^2+125}{a^2+25},\\ Det=\frac{a^2{d}^2\epsilon k^4+5abd\epsilon k^2+25d^2\epsilon }{a^2+25}\\ +\frac{k^4-4a^{2}d{k}^2+25ab+100d{k}^2+d{k}^2{a}^2+25d{k}^2}{a^2+25}.\end{array}$ (2.3)

As often as simple, we discuss the turing instability in the sequel paper which satisfy with the following two assumption,

H₀: The parameters b lying in the regime above the curve $5ab-3a^2+125>0$ ;

H₁: The necessary condition $-15625-9a^4+130a^{3}b\epsilon +\left(-25b^2{\epsilon }^2+750\right)a^2+1250ab\epsilon <0$ is satisfied since $De{t}_{\min }\left(d\right)<0$ .

To further compute the turing point, we solve $Det=0$ to get turing curves

$\begin{array}{l}{d}_{k\mathrm{,0}}=\frac{1}{2}\frac{-5ab\epsilon +3a^2-125+\sqrt{25a^2{b}^2{\epsilon }^2-130a^{3}b\epsilon +9a^4-1250ab\epsilon -750a^2+15625}}{\epsilon \left(a^2+25\right)k^2{\pi }^2}\mathrm{,}\\ d_{k\mathrm{,1}}=\frac{1}{2}\frac{-5ab\epsilon +3a^2-125+\sqrt{25a^2{b}^2{\epsilon }^2-130a^{3}b\epsilon +9a^4-1250ab\epsilon -750a^2+15625}}{\epsilon \left(a^2+25\right)k^2{\pi }^2}\end{array}$ (2.4)

The turing-turing bifurcation also occurs at $\left(b^{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{d}^{\mathrm{<mo>\; *\; </mo>}}\right)$ by setting $d_{k+1,0}=d_{k,1}$ to have

$\begin{array}{l}{b}_k^{*}=\frac{5\left(2a^2{k}^2+2a^{2}k+a^2+50k^2\right)\left(2k^2+2k+1\right)}{10k^2{\left(k+1\right)}^{2}a\epsilon }+\frac{3a^2-125}{5a\epsilon }\\ +\frac{2\sqrt{10}\sqrt{\left(a^2+25\right)\left(k^4+2k^3+\frac{13}{8}{k}^2+\frac{5}{8}k+\frac{5}{32}{a}^2+\frac{125}{8}{\left(k+\frac{1}{2}\right)}^2\right)}\left(2k^2+2k+1\right)}{5k^2{\left(k+1\right)}^{2}a\epsilon }\end{array}$ (2.5)

By the above discussion, $Det=0$ has a simple threshold curve $d=d\left(\epsilon \right)$ of Turing bifurcation as depicted by Equation (2.4) if and only if $b_{\ast ,k}<b<b_k^{\ast }$ . We also compute the transversal condition

$\frac{\text{d}\lambda }{\text{d}b}=\frac{-k^2{\pi }^2\left(a^2+25\right)\left(\epsilon +1\right)d\left(\epsilon \right)+3a^2-5ab-125}{k^2{\pi }^2\left(\left(a^2+25\right)d\left(\epsilon \right)k^2{\pi }^2+5ab\right)d\left(\epsilon \right)}<0$ (2.6)

Therefore we obtain the following results:

Proposition 2.1 With fixed parameter $a\mathrm{,}b\mathrm{,}\epsilon$ , the sequel points tracking on the Turing bifurcation curves has $d_{k,0}<d_{k,1}$ by Equation (2.4), for all positive integer k. However, turing-turing bifurcation occurs if and only if $d_{k+1,0}=d_{k,1}$ for some $k>1$ , therefore, the stability property of the homogeneous solution of system (1.1) (or system (2.1)), is changed when parameters pass over the Turing bifurcation curves, and the Turing instability regions are partitioned from the stability region.

For example, choosing $a=20,b=12$ , and by satisfying the basic assumptions H₀ and H₁, the turing bifurcation curves are drawn as shown in Figure 1(a) whilst $\epsilon$ lying below ${\epsilon }^{*}=0.91365064e-1$ . It is seen that the steady state is usually asymptotically stable if $\epsilon >{\epsilon }^{*}$ . The alike conclusion is satisfied for $a=20,\epsilon =0.05565$ , and we get the turing instability region given that the parameter b below the line $b=b^{*}=19.70136134$ . The Turing bifurcation curves are plotted as shown in Figure 1(b).

3. Pitchfork Bifurcation Branch from Turing Point

As discussed in section 2, if Turing condition is satisfied, the non-constant

Figure 1. The turing bifurcation of parabolic Equation (2.1). (a) The Turing bifurcation curve on d-ε plane with fixed parameters as $a=20,b=12$ . (b) The Turing bifurcation curve on b-d plane with fixed parameters as $a=20,\epsilon =0.05565$ . The interaction point which satisfies Equations (2.4) and (2.5) is happened with turing-turing bifurcation of codimension 2 sigularity.

steady solution arise from the Turing point. We give a description of Pitchfork solution branch bifurcating from Turing point which is also verified by germ’s strong equivalent property. To clarify the bifurcating solution classifying problem, we conclude the following proposition:

Proposition 3.1 Suppose on some neighborhood ${\mathcal{B}}_{x\mathrm{,}\lambda }$ of the trivial solution, we define a function $g\mathrm{<mo>\; :\; </mo>}{R}^2\to R$ which is $C^{\infty }$ , then a germ $g\in {\mathcal{B}}_{x\mathrm{,}\lambda }$ is strongly equivalent to polynomial $\alpha x^3\pm \beta \lambda x$ if and only if at $\left(x\mathrm{,}\lambda \right)=\left(\mathrm{0,0}\right)$ ,

$\begin{array}{l}g=\frac{\partial g}{\partial x}=\frac{{\partial }^{2}g}{\partial x^2}=\frac{\partial g}{\partial \lambda }=0\\ \alpha =sign\frac{{\partial }^{3}g}{\partial x^3},\beta =sign\frac{{\partial }^{2}g}{\partial \lambda \partial x}\end{array}$

Suppose $\alpha ={\alpha }_0$ is the bifurcation point of singularity. To discuss the stability property of the homogeneous steady state $\left(u_{\mathrm{<mo>\; *\; </mo>}}\mathrm{,}{v}_{\mathrm{<mo>\; *\; </mo>}}\right)$ , by doing axis transformation $u=u_1+u_{*}$ , $v=u_2+v_{*}$ , $\alpha ={\alpha }_0+\lambda$ , the corresponding parabolic system is written as the addition of the linear system $L\left(u_1\mathrm{,}{v}_1\mathrm{,}\lambda \right)$ and the corresponding nonlinear term $N\left(u_1\mathrm{,}{v}_1\mathrm{,}\lambda \right)$ , that is

$G\left(u_1\mathrm{,}{v}_1\mathrm{,}\lambda \right)=L\left(u_1\mathrm{,}{v}_1\mathrm{,}\lambda \right)+N\left(u_1\mathrm{,}{v}_1\mathrm{,}\lambda \right)$ (3.1)

with

$L\left(u_1,v_1,\lambda \right)=\left(\begin{array}{c}\Delta u_1\\ \Delta v_1\end{array}\right)+\left(\begin{array}{cc}\frac{3a^2-125}{a^2+25}& \frac{-20a}{a^2+25}\\ \frac{2b{a}^2}{a^2+25}& \frac{-5ba}{a^2+25}\end{array}\right)\left(\begin{array}{c}{u}_1\\ v_1\end{array}\right)$ (3.2)

and

$\begin{array}{l}N\left(u_1\mathrm{,}{v}_1\mathrm{,}\lambda \right)\\ =\left(\begin{array}{r}\hfill -\frac{20a\left(a^2-75\right)}{{\left(a^2+25\right)}^2}{u}_1^2+\frac{100a^2-2500}{{\left(a^2+25\right)}^2}{u}_1{u}_2+\frac{100a^4-15000a^2+62500}{{\left(a^2+25\right)}^3}{u}_1^3-\frac{500a\left(a^2-75\right)}{{\left(a^2+25\right)}^3}{u}_1^2{u}_2+o{\left(u_1+u_2\right)}^3\\ \hfill -\frac{5ba\left(a^2-75\right)}{{\left(a^2+25\right)}^2}{u}_1^2+\frac{25a^{2}b-625b}{{\left(a^2+25\right)}^2}{u}_1{u}_2+\frac{25b\left(a^4-150a^2+625\right)}{{\left(a^2+25\right)}^3}{u}_1^3-\frac{125ba\left(a^2-75\right)}{{\left(a^2+25\right)}^3}{u}_1^2{u}_2+o{\left(u_1+u_2\right)}^3\end{array}\right)\end{array}$ (3.3)

To compute the norm form near $\left(\mathrm{0,0}\right)$ , considering the linear operator $L\left(u_1\mathrm{,}{v}_1\mathrm{,}\lambda \right)$ , the corresponding Fredholm operator should has index zero, and we have the following proposition:

Proposition 3.2 Underlying the sigularity condimension 1 bifurcation at Turing point $\alpha ={\alpha }_0$ in system (1.2), the normal form of Equation (2.7) undertakes its strong equivalent form which can be expressed as $g=x^3+\lambda x$ . In addition, the Turing bifurcation is a supercitial Pitchfork if $\lambda <0$ , or either manifies a subcritical Pitchfork if $\lambda >0$ .

Hence after the recognition of normal form problem arising from Turing bifurcation with simple zero characteristic root is solved.

Proof: To carry out the normal form computation, we split space as

$X=Ker\left(L\right)+M,Y=R\left(L\right)+N$ (3.4)

With $Ker\left(L\right)$ and $R\left(L\right)$ are respectively the kernal space and the range space of linear operator L.

For example $Ker\left(L\right)=\left\{\beta \varphi \mathrm{<mo>\; |\; </mo>}\beta \ne 0\right\}$ , $\varphi$ is the unit eigenvector with

$\varphi =\sqrt{2}\cos \left(k\pi x\right){\left(\begin{array}{c}d{k}^2{\pi }^2\left(a^2+25\right)+5ba\\ 2b{a}^2\end{array}\right)}^{\text{T}}$

On the center manifold, define the projection mapping $E\mathrm{<mo>\; :\; </mo>}Y\to N$ which can separate Equation (2.7) into

$EG\left(u_1\mathrm{,}{v}_1\mathrm{,}\lambda \right)=0$ (3.5)

and the corresponding bifurcation equation

$\left(I-E\right)G\left(u_1\mathrm{,}{v}_1\mathrm{,}\lambda \right)=0$ (3.6)

By the space decomposition, we write $\left(\begin{array}{c}{u}_1\\ v_1\end{array}\right)=x\varphi +W$ , substituting it into Equation (2.9) and noticed that $EL=0$ to get $g\mathrm{<mo>\; :\; </mo>}Ker\left(L\right)\times M\times R\to N$ ,

$\begin{array}{c}g\left(x\mathrm{,}\lambda \right)=\langle \psi \mathrm{,}G\left(x\varphi +W\right)\rangle =\langle \psi \mathrm{,}L\left(x\varphi +W\right)+N\left(x\varphi +W\right)\rangle \\ =\langle \psi \mathrm{,}{L}_{\alpha }{}_0\varphi \rangle x+\langle \psi \mathrm{,}N\left(x\varphi +W\right)\rangle \end{array}$ (3.7)

with definition of $L_{\alpha }={\left(L_{\epsilon }\mathrm{,}{L}_d\right)}^{\text{T}}$ while regard the bifurcation point as

$L_{\alpha }={\left(L_{\epsilon }\mathrm{,}{L}_d\right)}^{\text{T}}$

or

$L_{\alpha }={\left(L_{\epsilon }\mathrm{,}{L}_b\right)}^{\text{T}}$

alike the definition at the turing point $L_{\alpha }={\left(L_{\epsilon }\mathrm{,}{L}_b\right)}^{\text{T}}$ respectively.

It is easily to compute that

$L_b=\left(\begin{array}{cc}0& 0\\ \frac{2a^2}{a^2+25}& \frac{-5a}{a^2+25}\end{array}\right)$

or

$L_d=-k_0^2{\pi }^2\left(\begin{array}{cc}\epsilon & 0\\ 0& 1\end{array}\right)$

We write the bifurcation equation as

$h\left(x,\lambda \right)=\left(I-E\right)G\left(x\varphi +W\right)=LW+N\left(x\varphi +W\right)-\varphi \langle \psi ,N\left(x\varphi +W\right)\rangle =0$ (3.8)

Since the operator L is Fredholm index 0, $L\mathrm{<mo>\; :\; </mo>}M\to R\left(L\right)$ is invertible mapping. Hence by the implicit function theorem, Equation (3.8) determines the unique expression $W=W\left(x\mathrm{,}\lambda \right)$ , then by Equation (3.7) we get the reduction equation on the center manifold

$g\left(x\mathrm{,}\lambda \right)=\langle \psi \mathrm{,}N\left(x\varphi +W\left(x\mathrm{,}\lambda \right)\right)\rangle$ (3.9)

To verify the strong equivalent form of germ g in proposition 3.1, suppose $L^{\mathrm{<mo>\; *\; </mo>}}$ being adjoint operator of L, we can choose

$\psi =\frac{\sqrt{2}}{{\left(d{k}^2{\pi }^2\left(a^2+25\right)+5ba\right)}^2-40b{a}^3}\cos \left(k\pi x\right){\left(\begin{array}{c}d{k}^2{\pi }^2\left(a^2+25\right)+5ba\\ -20a\end{array}\right)}^{\text{T}}$

Underlying Neumamn boundary condition, we compute $g_{xx}=0$ , and

$\begin{array}{l}{g}_{x\lambda }=\langle \psi ,L_{\alpha }\varphi \rangle x\lambda ,\\ g_{xxx}=\langle psi,d^{3}N\left(0,0\right)\left(\varphi ,\varphi ,\varphi \right)\rangle x^3,\end{array}$ (3.10)

Therefore, the normal form in Proposition 3.2 is verified and the pitchfork bifurcation direction is determined by the sign of $\lambda =\frac{g_{x\lambda }}{g_{xxx}}$ .

4. Numerical Simulation

Based on the results in section 2 and section 3, we can compute the corresponding Turing point via varying free parameters. For example, by simple calculation, the Turing bifurcation happens at the homogeneous equilibrium solution with chosen parameters $a=20$ , $b=16$ , $\epsilon =0.05565$ whilst $k=1$ , $d=1.1181$ , and from the above discussion in section 3, we can compute the base $\varphi$ of kernal

of linear operator L as $\varphi =\cos \left(k\pi x\right)\left(\begin{array}{c}0.44103\\ 0.89749\end{array}\right)$ , which further derive the formula

$g\left(x\mathrm{,}{d}_{\epsilon }\right)=0.4001062699e-1x^3+0.8362161555d_{\epsilon }x$ , then the nonhomogeneous steady state solution branch bifurcates from Turing point in accordance with the direction $d_{\epsilon }<0$ , which is pitchfork and subcritical bifurcation. However with $k=2$ , $d=0.7764$ , we have $g\left(x,d_{\epsilon }\right)=0.1050273688x^3-1.282881210d_{\epsilon }x$ , which is supercritical pitchfork of Turing bifurcation. Choosing $a=20$ , $\epsilon =0.05565$ to get Turing point $k=1$ , $d=1$ , $b=14.96776630$ and $k=2$ , $d=0.83$ , $b=14.40714340$ , the normal form at two different bifurcation point are respectively $g\left(x,b_{\epsilon }\right)=0.3884830554e-1x^3-0.1144224171b_{\epsilon }x$ and $g\left(x,b_{\epsilon }\right)=0.1191247571x^3-0.4528851324e-1b_{\epsilon }x$ , which is Turing subcritical bifurcation. As shown in Figure 2(a), with fixed parameters $a=20$ , $\epsilon =0.05565$ , $b=16$ , the continuation of Turing bifurcation is carrying out with the continuous varying of parameter d, which illustrates the Turing bifurcation point at $d=1.1181$ which is subcritical. As shown in Figure 2(b), the parameters are fixed as $a=20$ , $\epsilon =0.05565$ , $d=1$ , we vary parameter b continuously, the continuation of supercritical Turing bifurcation at $b=14.96776630$ is simulated. The simulation algorithms are as often familiar as the well known differential quadrature algorithms, element free Galerkin method, and trigonometric B-spline functions interpolation method, etc.

The temporal-spatial solutions near Turing points as shown in Figure 3 illustrate the direction of pitchfork bifurcation, which is in coincidence with the sign of the coefficients $\lambda$ in normal form. For example, fixed parameter with $a=20$ , $b=16$ , $\epsilon =0.05565$ , as shown in Figure 3(a) and Figure 3(b), the subcritical Turing bifurcation happens at $d=1.1181$ . The constant steady state is asymptotically stable at $d=1.1190$ and the non-homogeneous solution is observed at $d=1.1080$ . The supercritical Turing bifurcation manifests bifurcating non-constant solution at $d=0.7804$ as shown in Figure 3(c), however the constant steady state is stable at $d=0.7704$ as shown in Figure 3(d). The Turing solution and the Turing oscillation solution are observed at $\tau =0.05$ , $d=0.78$ and $\tau =0.07$ , $d=1.112$ , as shown in Figure 3(e) and Figure 3(f), respectively.

Figure 2. The continuation of turing bifurcation with $a=20$ , $\epsilon =0.05565$ . (a) Chosen with $b=16$ , the subcritical turing bifurcation happens at $d=1.1181$ ; (b) Chosen with $d=1$ , the supercritical turing bifurcation happens at $b=14.96776630$ .

Figure 3. The supercritical turing bifurcation with $a=20$ , $b=16$ , $\epsilon =0.05565$ and (a) $d=1.1190$ , (b) $d=1.1080$ ; The subcritical turing bifurcation with (c) $d=0.7804$ ; (d) $d=0.7704$ ; (e) $\tau =0.05$ , $d=0.78$ (f) $\tau =0.07$ , $d=1.112$ .

Figure 4. The turing bifurcation and turing oscillation occurs via varying time delay, the parameters are fixed with $a=20$ , $b=12$ , $\epsilon =0.04565$ , $d=0.89$ . (a) Nonhomogeneous solution occurs with $\tau =0$ ; (b) The turing oscillation solution with $\tau =0.04$ ; (c) The turing oscillation solution with $\tau =0.049$ ; (d) The turing oscillation solution with $\tau =0.0498$ ; (e) The wave pattern with $\tau =0.049$ ; (f) The wave pattern with $\tau =0.0498$ .

We also simulate the 2D tempo-spatial solutions and produce pictures via stretching y-axis in accordance with x-axis direction, which show us the numerical solutions of PDE in an easy way, as shown in Figures 4(a)-(d). With fixed parameter $a=20$ , $b=12$ , $\epsilon =0.04565$ , $d=0.89$ , varying time delay, system (1.1) manifests the nonhomogeneous solution at $\tau =0$ , however, the occurrence of the Turing oscillation are simulated at $\tau =0.04$ , $\tau =0.049$ and $\tau =0.0498$ , respectively. The periodical wave patterns are schemed by projection onto X-Y plane hence is given in Figure 4(e) and Figure 4(f).

We further do simulation to verify the wave pattern formation in Gray-Scott diffusion model with time delay effects. We choose $\tau =0.4$ , fixed parameters with $a=20$ , $b=9.6$ , the wave pattern formation is simulated which illustrates the route of periodical oscillation to chaos, as shown in Figures 5(a)-(d), With different value of $\epsilon$ and d, from its periodicity to irregularity, the wave patterns underlying Turing oscillation are simulated, which manifests turing-turing bifurcation can bring complex dynamical behavior underlying time delay effects.

Figure 5. The wave pattern formation with $a=20$ , $b=9.6$ and (a) $\epsilon =0.02$ , $d=0.4565$ , (b) $\epsilon =0.018$ , $d=0.3565$ ; (c) $\epsilon =0.015$ , $d=0.35$ ; (d) $\epsilon =0.01565$ , $d=0.3$ .

5. Discussion

We discuss the dynamics in Gray-Scott diffusion model produced via turing bifurcation, which usually emphasis on the significant turing condition. In an obvious way, the turing-turing bifurcation can be either independent on time delay. Underlying small time delay, the Gray-Scott diffusion model brings forth the rich formation of wave patterns. By applying Lyapunov-Schimdt reduction method, the normal form of the turing bifurcation was computed. Correspondingly, the pitchfork bifurcation direction is determined by the coefficients of the strong equivalent form of germs in R². However, the turing oscillation was observed as ascending time delay, which discovers that system of PDE may manifest codimesion-2 singularity and the bifurcation mechanism need to be further studied.

Acknowledgements

Sincere thanks to the members of IJMNTA for their professional performance, and special thanks to managing editor for a rare attitude of high quality.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Smith, H. (2011) An Introduction to Delay Differential Equations with Applications to the Life Sciences. Springer, Berlin. https://doi.org/10.1007/978-1-4419-7646-8
[2]	Kuang, K. (1992) Delay Differential Equations with Application in Population Dynamics. Springer, New York.
[3]	Skalski, G. and Gilliam, J.F. (2002) Functional Responses with Predator Interference: Viable Alternatives to the Holling Type II Model. Ecology, 82, 3083-3092. https://doi.org/10.1890/0012-9658(2001)082[3083:FRWPIV]2.0.CO;2
[4]	Jeschke, J., Kopp, J.M. and Tollrian, R. (2002) Predator Functional Responses: Discriminating between Handling and Digesting Prey. Ecological Monographs, 72, 95-112. https://doi.org/10.1890/0012-9615(2002)072[0095:PFRDBH]2.0.CO;2
[5]	Chang, S.L. and Jeng, B.W. (2007) Liapunov Schmidt Reduction and Continuation for Nonlinear Schordinger Equations. SIAM Journal on Scientific Computing, 29, 729-755. https://doi.org/10.1137/050642861
[6]	Bohmer, K. (1993) On a Numerical Liapunov-Schmidt Method for Operator Equations. Computing, 51, 237-269. https://doi.org/10.1007/BF02238535
[7]	Bohmer, K. and Mei, Z. (1990) On a Numerical Liapunov-Schmidt Method. In: Allgower, E.L. and Georg, K., Eds., Computational Solution of Nonlinear Systems of Equations, Lectures in Applied Mathematics, Vol. 26, 79-98.
[8]	Ashwin, P., Bohmer, K. and Mei, Z. (1995) A Numerical Liapunov-Schmidt Method with Applications to Hopf Bifurcation on a Square. Mathematics of Computation, 64, 649-670. https://doi.org/10.1090/S0025-5718-1995-1284661-X
[9]	Shi, H.B., Li, W.T. and Lin, G. (2010) Positive Steady States of a Diffusive Predator Prey System with Modified Holling Tanner Functional Response. Nonlinear Analysis: Real World Applications, 11, 3711-3721. https://doi.org/10.1016/j.nonrwa.2010.02.001
[10]	Li, M.F., Han, B., Xu, L. and Zhang, G. (2013) Spiral Patterns Near Turing Instability in a Discrete Reaction Diffusion System. Chaos, Solitons & Fractals, 49, 1-6. https://doi.org/10.1016/j.chaos.2013.01.010
[11]	Mukhopadhyay, B. and Bhattacharyya, R. (2006) Modeling the Role of Diffusion Coefficients on Turing Instability in a Reaction Diffusion Prey Predator System. Bulletin of Mathematical Biology, 68, 293-313. https://doi.org/10.1007/s11538-005-9007-2
[12]	Hsu, S.B. and Huang, T.W. (1995) Global Stability for a Class of Predator Prey Systems. SIAM Journal on Applied Mathematics, 55, 763-783. https://doi.org/10.1137/S0036139993253201
[13]	Fedoseyev, A.I., Friedman, M.J. and Kansa, E.J. (2002) Improved Multiquadric Method for Elliptic Partial Differential Equations via PDE Collocation on the Boundary. Computers & Mathematics with Applications, 43, 439-455. https://doi.org/10.1016/S0898-1221(01)00297-8
[14]	Hioe, F.T. (1999) Solitary Waves for N Coupled Nonlinear Schrodinger Equations. Physical Review Letters, 82, 1152-1155. https://doi.org/10.1103/PhysRevLett.82.1152
[15]	Ahmed, B. (2012) 1-Soliton Solution for Zakharov Equations. Advanced Studies in Theoretical Physics, 6, 457-466.
[16]	Bao, W. and Yang, L. (2007) Efficient and Accurate Numerical Methods for the Klein-Gordon-Schrödinger Equations. Journal of Computational Physics, 225, 1863-1893. https://doi.org/10.1016/j.jcp.2007.02.018
[17]	Hu, G.P., Li, X.L. and Wang, Y.P. (2015) Pattern Formation and Spatiotemporal Chaos in a Reaction—Diffusion Predator—Prey System. Nonlinear Dynamics, 81, 265-275. https://doi.org/10.1007/s11071-015-1988-2
[18]	Banerjee, M. and Banerjee, S. (2012) Turing Instabilities and Spatio-Temporal Chaos in Ratio-Dependent Holling-Tanner Model. Mathematical Biosciences, 236, 64-76. https://doi.org/10.1016/j.mbs.2011.12.005