Turing instability in a modified cross-diffusion Leslie–Gower predator–prey model with Beddington–DeAngelis functional response

In this paper, a modified cross-diffusion Leslie–Gower predator–prey model with the Beddington–DeAngelis functional response is studied. We use the linear stability analysis on constant steady states to obtain sufficient conditions for the occurrence of Turing instability and Hopf bifurcation. We show that the Turing instability and associated patterns are induced by the variation of parameters in the cross-diffusion term. Some numerical simulations are given to illustrate our results.


Introduction
In this paper, we consider the modified Leslie-Gower predator-prey model where is a bounded domain in R N with smooth boundary ∂ , u and v represent the densities of the prey and the predator, is the Laplacian operator, and ν is the outward unit normal vector on ∂ . Also r is the growth rate of the prey in the absence of the predator, k is the carrying capacity of the prey in the absence of the predator population, d is the intrinsic growth rate of the predator species, e is the maximum rate of death of the predator population, and μ > 0 is the linear diffusion coefficient. The nonlinear cross diffusion term α (uv) describes a tendency such that the prey species keep away from high-density areas of the predator species and ( v 1+βu ) models a situation in which the population pressure of the predator species weakens in the high density place of the prey species. This model is called the modified Leslie-Gower predator-prey model. For more information about the Leslie-Gower predator-prey model, we refer the reader to [9][10][11]. The term p(u, v) = au m + bu + cv is called the Beddington-DeAngelis functional response. Here a, b, c, and m are the consumption rate, the saturation constant for an alternative prey, the predator interference, and the coefficient of environmental protection for the prey, respectively. To find more details on the background of this functional response, see [1,4]. In this paper we assume that all constants a, b, c, d, e, k, m, r, α, and β are positive.
As far as we know, there are few contributions on the cross-diffusion prey-predator model with the nonlinear cross diffusion term ( v 1+βu ). For example, in [3] the authors considered the following predator-prey model: x ∈ , t > 0, ∂u 2 ∂t = ((μ + 1 1+βu 1 )u 2 ) + u 2 (bu 2 + du 1 1+mu 1 ), x ∈ , t > 0, u 1 (x, 0) = u 2 (x, 0) ≥ 0, x ∈ , They investigated the existence of the steady-state bifurcation and the Turing instability using the Lyapunov-Schmidt reduction. In [22], the authors proved the existence and multiplicity of positive steady state solutions for system (1.2). They also established the uniqueness of positive solution in a special domain. Yan and Guo in [20] discussed the direction and stability of the Hopf bifurcation for a Lotka-Volterra model with delay and cross-diffusion by the Lyapunov-Schmidt reduction method. Further, the effect of delay on the stability and Turing instability of equilibrium points for the Crowley-Martin predator-prey model was studied in [6]. In [12], the existence of a positive periodic solution for a nonautonomous Leslie-Gower reaction-diffusion model has been proved. The authors used the comparison theory of differential equations to get their results. Surendar et al. [17] studied the local stability of positive equilibrium and the existence of a Hopf bifurcation for a diffusive Holling-Tanner predator-prey system with the stoichiometric density dependence. The authors in [2] derived the Hopf bifurcation and its properties for a class of system of reaction-diffusion equations with two discrete time delays. In [21], Yi et al. studied the Turing instability and the Hopf bifurcation for the Lengyel-Epstein system. Many research works have been devoted to the Hopf bifurcation for ODE preypredator models. For instance, Rihan et al. studied the stability of the Hopf bifurcation for a fractional-order delay prey-predator system in [15]. Sivasamy et al. in [16] analyzed the direction and stability of a Hopf bifurcation for a modified Leslie-Gower prey-predator model with gestation delay. In [13] the authors discussed the local and global stability of the interior equilibrium point and the Hopf bifurcation for a prey-predator model with Hassell-Varley functional response.
The main purpose of this paper is to investigate the effect of the cross-diffusion term on the stability and pattern formation in system (1.1). The existence of a periodic solution in the instability region and the Hopf bifurcation for the ODE and PDE models are also studied. We discuss the stability and direction of periodic solutions for the ODE system by the Poincare-Andronov-Hopf bifurcation theorem. For the PDE system, we use the center manifold theorem and the normal form theory. Our numerical results show that a Turing pattern is emerged by the variation of cross-diffusion parameters.
The organization of the rest of the paper is as follows: In Sect. 2, we analyze the asymptotic behavior of stationary solutions and the Hopf bifurcation for the ODE system. In Sect. 3, the Turing instability, the stability of bifurcating periodic solutions, and the bifurcation direction are investigated for system (1.1). Some numerical simulations are presented in Sect. 4.

Stability analysis and Hopf bifurcation for the ODE system
The ODE system corresponding to system (1.1) is as follows: We consider d as the bifurcation parameter. System (2.1) has the following equilibrium points: , (m + bk)r(e + c)abk er(e + c) .
We can rewrite v * as v * = m e + bk e 1 -a r(e + c) .
Hence a sufficient condition for the positivity of u * and v * is From now on we assume that (2.2) is satisfied. The Jacobian matrix of system (2.1) at U * is and the characteristic equation of J is given by where T 0 = trac(J(U * )) and D 0 = det(J(U * )). Since (u * , v * ) is an equilibrium of system (2.1), we can simplify T 0 and D 0 as follows: Then the roots of (2.4) are λ = α 0 (d) ± iβ 0 (d) where The root of the equation T 0 (d) = 0 is Under the condition d 0 is positive. Therefore, by the fact that D 0 > 0, we conclude that J(U * ) has a pair of imaginary eigenvalues at d = d 0 as Also, Then, by the Hopf theorem [8], system (2.1) has a Hopf bifurcation at U * when d = d 0 . Furthermore, T 0 = d 0d, then T 0 > 0 for 0 < d < d 0 , which implies that the characteristic equation (2.4) has a root with positive real part. So the equilibrium solution U * is unstable. Also T 0 < 0 for d 0 < d, therefore both roots of (2.4) have negative real parts. Hence, the equilibrium solution U * is locally asymptotically stable.
In the sequel, we detail the direction of Hopf bifurcation and the stability of bifurcating periodic solutions for system (2.1). First, we translate U * to the origin. After this translation, we obtain the following system: where where is a multi-index with | | = 4. Set Consider the following transformation: So the Jordan canonical form of system (2.7) is as follows: and After transforming (2.10) into the normal form and rewriting it in the polar coordinates, we obtain the following equations: For more details about the normal form in polar coordinates (2.11), see Sect. 19.2 of [19].
Using Taylor expansion of the coefficients in (2.11) about d = d 0 , we geṫ By Lemma 20.2.2 in [19], the sign of a(d 0 ) determines the stability of the periodic orbit. The coefficient a(d 0 ) is given by By calculation, we have Then we get (2.12) (2.13) (2.14) Then a(d 0 ) = X 1 + X 2 + X 3 16 .
If (bk + m) r(e + c)a 2ma 2 > 0, (2.15) then by (2.2) and (2.5), X 1 , X 2 , and X 3 are negative. So a(d 0 ) < 0. By the above calculation and the Poincare-Andronov-Hopf bifurcation theorem [19], we obtain the following result. According to Theorem 2.1, system (2.1) has periodic solutions bifurcating from the Hopf bifurcation. In the next theorem, we show that system (2.1) has at least one stable periodic solution for d < d 0 .
Theorem 2.2 Let a < rc and (2.5) be satisfied. For d < d 0 , system (2.1) has at least one stable periodic solution U(t) = (u(t), v(t)) such that ε < u(t) < k and m e < v(t) < m+bk e for 0 < ε < 1 -a rc k. (2.16) Proof We construct a trapping region satisfying the conditions of the Poincare-Bendixson theorem. To do this, let The segments C 1 , C 2 , C 3 , and C 4 form a Jordan curve C. Denote the interior of C by D. One can see U * = (u * , v * ) is the only positive equilibrium point in D. On the other hand, dv dt | (u,v)∈C 1 > 0, du dt | (u,v)∈C 2 < 0, and dv dt | (u,v)∈C 3 < 0. Also, by (2.16), for (u, v) ∈ C 4 , we have .
Hence du dt | (u,v)∈C 4 > 0. Therefore on the boundary of D the vector filed points inwards the set. This implies that a trajectory will stay in D for all t once it has entered the set. Since U * = (u * , v * ) ∈ D is unstable, by the Poincare-Bendixson theorem, system (2.1) has at least one stable periodic solution in D.

Turing instability of the PDE system
In this section, we obtain a condition under which the cross-diffusion system (1.1) loses stability and the stable positive equilibrium becomes Turing unstable.
Let {μ n } ∞ n=0 be the eigenvalues of -in with zero Neumann boundary conditions. The linearized system of (1.1) at U where J is defined by (2.3) and Then λ is an eigenvalue of L := D + J(U * ) on if and only if λ is an eigenvalue of the matrix -μ n D + J(U * ) for n ≥ 0. Therefore λ satisfies the characteristic equation where According to (3.4), In the next theorem, we summarized the asymptotic behavior of the equilibrium U * for system (1.1).
(ii). Since ev * m+bu * = 1, from (1) we conclude Furthermore, D 0 (d) > 0 for d > 0, and we can rewrite D n (d) as By (3.7), we have D n (d) > 0 for n ∈ N and d > 0. Since T n (d) < 0 for n ∈ N ∪ {0} when d > d 0 , we deduce that the real part of every eigenvalue of L is negative. Therefore, U * is locally asymptotically stable for system (1.1). Now, let (2) be satisfied. Then D 1 = αv * N 1 + M 1 > 0. Also, we can rewrite M n as Hence, for n = 2, 3, . . . , using condition (3.6) we have Since N n+1 ≥ N n for n ∈ N, the above inequalities imply D n (d) ≥ D 1 (d) > 0 for n ≥ 2 and d > d 0 . Again using T n (d) < 0 for n ∈ N, we conclude that all roots of the characteristic equation (3.2) have negative real part for n ≥ 0 when d > d 0 . Therefore U * is locally asymptotically stable for system (1.1).
The linearized system of (3.9) at U * has the following form: where J and D are defined by (2.3) and (3.1) respectively and Similarly, by transferring the equilibrium point U * to the origin, system (3.9) becomes where f 1 and f 2 are defined by (2.8) and (2.9). Define the operator L * by where J * is the conjugate of J. Then L * is the adjoint operator of L. Set where β 0 (d 0 ) = ru * d 0 k and i is the imaginary unit. Then T g dx is the inner product in L 2 (0, lπ) × L 2 (0, lπ). We decompose X as According to [8], for every U T = (u, v) ∈ X, we have where zq +zq ∈ X c and w T = (w 1 , w 2 ) ∈ X s . Then Thus system (3.12) in (z, w) coordinates is as follows: H(z,z, w), (3.14) where f ∨ = (f 1 , f 2 ) T is defined by (2.8)-(2.9) and By calculation, we obtain Then H(z,z, w) = 0. Besides, the center manifold of system (3.14) is illustrated as follows: Then we get So w 20 = w 02 = w 11 = 0. Therefore, the restriction of (3.14) to the center manifold in z and z coordinates is given by where h 20 = q * , B(q, q)  can see the mentioned result in Fig. 1(b) where the initial conditions are taken at (0.03, 1.3) and (0.05, 1.4). In [7], we can see the description of periodic orbits by using the fractal-like behavior of prime numbers (also, see [14,18]). Now, consider system (1.1) in = (0, 16π) × (0, 16π) with parameters defined by (4.1) and μ = 103, β = 87.   Also, we derive a(d 0 ) = -0.0002. By Theorem 2.1, the positive equilibrium (u * , v * ) is asymptotically stable when d > d 0 and unstable when d < d 0 . Moreover, a Hopf bifurcation occurs at d = d 0 . The direction of Hopf bifurcation is subcritical and the bifurcating periodic solutions are asymptotically orbitally stable. For d = 0.5 > d 0 , the equilibrium (u * , v * ) is asymptotically stable and the solution goes to (u * , v * ). This result is shown in Fig. 3(a) where the initial condition is taken at (0.4, 5.5). In Fig. 3(b), we choose d = 0.45 < d 0 , and