Inhomogeneous spatial patterns in diffusive predator-prey system with spatial memory and predator-taxis

In this paper, we consider a diffusive predator-prey system with spatial memory and predator-taxis. Since in this system, the memory delay appears in the diffusion term, and the diffusion term is nonlinear, the classical normal form of Hopf bifurcation in the reaction-diffusion system with delay can't be applied to this system. Thus, in this paper, we first derive an algorithm for calculating the normal form of Hopf bifurcation in this system. Then in order to illustrate the effectiveness of our newly developed algorithm, we consider the diffusive Holling-Tanner model with spatial memory and predator-taxis. The stability and Hopf bifurcation analysis of this model are investigated, and the direction and stability of Hopf bifurcation periodic solution are also researched by using our newly developed algorithm for calculating the normal form of Hopf bifurcation. At last, we carry out some numerical simulations, two stable spatially inhomogeneous periodic solutions corresponding to the mode-1 and mode-2 Hopf bifurcations are found, which verifies our theoretical analysis results.


Introduction
The reaction-diffusion equations based on the Fick's law have been widely used in physics, chemistry and biology [1,2,3].More precisely, based on the Fick's law, that is the movement flux is in the direction of negative gradient of the density distribution function, the diffusive Brusselator model with gene expression time delay [4], the diffusive predator-prey model in heterogeneous environment [5,6], the diffusive predatorprey model with a protection zone [7], the diffusive predator-prey model with prey social behavior [8], the diffusive predator-prey model with protection zone and predator harvesting [9] have been studied by many scholars.Furthermore, in order to include the episodic-like spatial memory of animals, Shi et al. [10] directed movement toward the negative gradient of the density distribution function at past time, and they proposed the following diffusive model with spatial memory where u(x, t) is the population density at spatial location x and time t, d 1 and d 2 are the Fickian diffusion coefficient and the memory-based diffusion coefficient, respectively, Ω ⊂ R is a smooth and bounded domain, u 0 (x, t) is the initial function, ∆u(x, t) = ∂ 2 u(x, t)/∂x 2 , u x (x, t) = ∂u(x, t)/∂x, u x (x, t−τ ) = ∂u(x, t−τ )/∂x, u xx (x, t − τ ) = ∂ 2 u(x, t − τ )/∂x 2 , and n is the outward unit normal vector at the smooth boundary ∂Ω.
Here, the time delay τ > 0 represents the averaged memory period, which is usually called as the memory delay, and f (u(x, t)) describes the chemical reaction or biological birth/death.Moreover, in order to further investigate the influence of memory delay on the stability of the positive constant steady state, on the basis Yehu Lv School of Mathematical Sciences, Beijing Normal University, Beijing 100875, China E-mail: mathyehul@163.com of model (1.1), Shi et al. [11] studied the spatial memory diffusion model with memory and maturation delays.Furthermore, Song et al. [12] proposed the following memory-based diffusion system subjects to homogeneous Neumann boundary condition ∂u(x, t) ∂t = d 11 ∆u(x, t) + f (u(x, t), v(x, t)) , x ∈ Ω, t > 0, ∂v(x, t) ∂t = d 22 ∆v(x, t) − d 21 (v(x, t)u x (x, t − τ )) x + g (u(x, t), v(x, t)) , x ∈ Ω, t > 0, u x (0, t) = u x (ℓπ, t) = v x (0, t) = v x (ℓπ, t) = 0, t ≥ 0, u(x, t) = u 0 (x, t), x ∈ Ω, − τ ≤ t ≤ 0, v(x, t) = v 0 (x), x ∈ Ω, where Ω = (0, ℓπ) with ℓ ∈ R + , u(x, t) and v(x, t) represent the densities of prey and predator at the location x and time t, respectively, d 11 ≥ 0 and d 22 ≥ 0 are the random diffusion coefficients, d 21 ≥ 0 is the memorybased diffusion coefficient, v 0 (x) is also the initial function, and f (u(x, t), v(x, t)) and g (u(x, t), v(x, t)) are the reaction terms.Furthermore, they derived an algorithm for calculating the normal form of Hopf bifurcation in the system (1.2).By numerical simulations, the stable spatially inhomogeneous periodic solutions and the transition from the unstable mode-2 spatially inhomogeneous periodic solution to the stable mode-1 spatially inhomogeneous periodic solution are found.
For the general predator-prey models in ecology, in addition to random diffusion of the predator and prey populations, the spatial movement of predator and prey populations also occurs, which is usually shown by the predator pursuing prey and prey escaping from predator [13].The pursuit and evasion between the predator and prey populations also have a strong impact on the movement pattern of the predator and prey populations [14,15,16].Furthermore, by noticing that such movement is not random but directed, i.e., predator moves toward the gradient direction of prey distribution, which is called prey-taxis, or prey moves opposite to the gradient of predator distribution, which is called predator-taxis.Recently, the predator-prey model with prey-taxis [17,18,19,20,21], the predator-prey model with indirect prey-taxis [22,23,24], the predator-prey model with predator-taxis [25] and the predator-prey model with indirect predator-taxis [26] have been researched.Especially, Wang et al. [13] considered the following diffusive predator-prey model with both predator-taxis and prey-taxis.The corresponding mathematical model is ∂u(x, t) ∂t = d∆u(x, t) + ξ (u(x, t)v x (x, t)) x + f (u(x, t), v(x, t)), x ∈ Ω, t > 0, where u 0 (x) is the initial function, d is the rescaled diffusion coefficient for the prey population, and the diffusion coefficient of the predator population is rescaled as 1.Furthermore, the term ξ (u(x, t)v x (x, t)) x represents the prey moves away from predator, and ξ ≥ 0 is the intrinsic predator-taxis rate.The term −η (v(x, t)u x (x, t)) x represents the predator moves towards prey, and η ≥ 0 is the intrinsic prey-taxis rate.Therefore, based on the models (1.2) and (1.3), and by considering that the spatial memory and predatortaxis, we proposed the following diffusive predator-prey model with spatial memory and predator-taxis subjects to homogeneous Neumann boundary condition x ∈ (0, ℓπ). (1.4) This paper is organized as follows.In Section 2, we derive an algorithm for calculating the normal form of Hopf bifurcation in the model (1.4).In Section 3, we obtain the normal form of Hopf bifurcation truncated to the third-order term by using our newly developed algorithm developed in Sec.2, and the mathematical expressions of the corresponding coefficients are given.In Section 4, we consider the diffusive Holling-Tanner model with spatial memory and predator-taxis.The stability and Hopf bifurcation analysis in this model are studied, and some numerical simulations are also given.In Section 5, we give a brief conclusion and discussion.
2 Algorithm for calculating the normal form of Hopf bifurcation in the system (1.4)

Characteristic equation at the positive constant steady state
Define the real-valued Sobolev space with the inner product defined by where the symbol T represents the transpose of vector, and let C := C([−1, 0]; X) be the Banach space of continuous mappings from [−1, 0] to X with the sup norm.It is well known that the eigenvalue problem where e j , j = 1, 2 is the unit coordinate vector of R 2 , and n ∈ N 0 = N ∪ {0} is often called wave number, N 0 is the set of all non-negative integers, N = {1, 2, ...} represents the set of all positive integers.Without loss of generality, we assume that E * (u * , v * ) is the positive constant steady state of system (1.4).The linearized equation of (1.4) where and Here, det(.)represents the determinant of a matrix, I 2 is the identity matrix of 2 × 2, and D 1 , D 2 , A are defined by (2.3).Then we can obtain where
Let τ = τ c + µ, |µ| ≪ 1 such that µ = 0 corresponds to the Hopf bifurcation value for system (1.4).Moreover, we shift E * (u * , v * ) to the origin by setting and normalize the delay by rescaling the time variable t → t/τ .Furthermore, we rewrite U (t) for U (x, t), and U t ∈ C for U t (θ) = U (x, t + θ), − 1 ≤ θ ≤ 0. Then the system (1.4) becomes the compact form where for ϕ = ϕ (1) , ϕ (2) T ∈ C, d(µ)∆ is given by x (0)ϕ x (0) + ϕ (1) (0)ϕ and F : C × R 2 → X is given by In what follows, we assume that F (ϕ, µ) is C k (k ≥ 3) function, which is smooth with respect to ϕ and µ.Notice that µ is the perturbation parameter and is treated as a variable in the calculation of normal form.Moreover, by (2.10), if we denote L 0 (ϕ) = τ c Aϕ(0), then (2.8) can be rewritten as where the linear and nonlinear terms are separated, and Thus, the linearized equation of (2.12) can be written as where Γ n (λ) = det M n (λ) with By comparing (2.16) with (2.5), we know that (2.15) has a pair of purely imaginary roots ±iω c for n = n c ∈ N, and all other eigenvalues have negative real parts, where ω c = τ c ω nc .In order to write (2.12) as an abstract ordinary differential equation in a Banach space, follows by [27], we can take the enlarged space then the equation (2.12) is equivalent to an abstract ordinary differential equation on BC Here, A is a operator from C 1 0 = {ϕ ∈ C : φ ∈ C, ϕ(0) ∈ dom(∆)} to BC, which is defined by and X 0 = X 0 (θ) is given by In the following, the method given in [27] is used to complete the decomposition of BC.
, where R 2 * is the two-dimensional space of row vectors, and define the adjoint bilinear form on C * × C as follows By choosing where the col(.)represents the column vector, φ(θ) = col (φ is the eigenvector of (2.14) associated with the eigenvalue iω c , and ψ(s) = col (ψ 1 (s), ψ 2 (s)) = ψe −iωcs ∈ C 2 with ψ = col (ψ 1 , ψ 2 ) is the corresponding adjoint eigenvector such that According to [27], the phase space C can be decomposed as where for φ ∈ C, the projection π : C → P is defined by (2.17) Therefore, by following the method in [27], BC can be divided into a direct sum of center subspace and its complementary space, that is where dim P = 2.It is easy to see that the projection π which is defined by (2.17), is extended to a continuous projection (which is still denoted by π), that is, π : BC → P. In particular, for α ∈ C, we have By combining with (2.18) and (2.19), U t (θ) can be decomposed as where w = col (w 1 , w 2 ) and If we assume that then (2.20) can be rewritten as Then by combining with (2.21), the system (2.12) is decomposed as a system of abstract ordinary differential equations (ODEs) on R 2 × Ker π, with finite and infinite dimensional variables are separated in the linear term.That is nc where I is the identity matrix, z = (z 1 , z 2 ) T , B = diag {iω c , −iω c } is the diagonal matrix, and

Consider the formal Taylor expansion
From (2.13), we have By combining with (2.19), the system (2.22) can be rewritten as where nc (2.25) In terms of the normal form theory of partial functional differential equations [27], after a recursive transformation of variables of the form where z, z ∈ R 2 , w, w ∈ Q 1 and U 1 j : R 3 → R 2 , U 2 j : R 3 → Q 1 are homogeneous polynomials of degree j in z and µ, a locally center manifold for (2.12) satisfies w = 0 and the flow on it is given by the two-dimensional ODEs which is the normal form as in the usual sense for ODEs.
By following [27] and [28], we have and where Proj p (q) represents the projection of q on p, and f 1 3 (z, 0, µ) is vector and its element is the cubic polynomial of (z, µ) after the variable transformation of (2.26), and it is determined by (2.38), and (2.30) In the following, for notational convenience, we let We then calculate g 1 j (z, 0, µ), j = 2, 3 step by step.
Furthermore, from (2.37), (2.39) and (2.40), we have and then we obtain Proj S (D w,wx,wxx f where 3 Normal form of the Hopf bifurcation and the corresponding coefficients According to the algorithm developed in Section 2, we can obtain the normal form of the Hopf bifurcation truncated to the third-order term ż = Bz + 1 2 where Here, B 1 is determined by (2.36), B 21 , B 22 and B 23 are determined by (2.42), (2.50), (2.51), respectively, and they can be calculated by using the MATLAB software.The normal form (3.1) can be written in real coordinates through the change of variables z 1 = v 1 − iv 2 , z 2 = v 1 + iv 2 , and then changing to polar coordinates by v 1 = ρ cos Θ, v 2 = ρ sin Θ, where Θ is the azimuthal angle.Therefore, by the above transformation and removing the azimuthal term Θ, (3.1) can be rewritten as where According to [29], the sign of K 1 K 2 determines the direction of the Hopf bifurcation, and the sign of K 2 determines the stability of the Hopf bifurcation periodic solution.More precisely, we have the following results (i) when K 1 K 2 < 0, the Hopf bifurcation is supercritical, and the Hopf bifurcation periodic solution is stable for K 2 < 0 and unstable for K 2 > 0; (ii) when K 1 K 2 > 0, the Hopf bifurcation is subcritical, and the Hopf bifurcation periodic solution is stable for K 2 < 0 and unstable for K 2 > 0.

Application to Holling-Tanner model with spatial memory and predator-taxis
In this section, we apply our newly developed algorithm in Section 2 to the Holling-Tanner model with spatial memory and predator-taxis, i.e., in the model (1.4), we let where β > 0, m > 0 and s > 0. Thus, (1.4) becomes the following model 2) The Holling-tanner model is one of the typical predator-prey models.For the ordinary differential equation (4.2) with d 11 = ξ = d 21 = d 22 = 0, it has been completely analyzed in [30].For the diffusive model (4.2) with ξ = d 21 = 0, the global stability of the positive constant steady state was proved in [31,32], and the Hopf bifurcation and Turing instability have been studied in [33].
In the following, we let Furthermore, let λ = iω n (ω n > 0) be a root of (4.5).By substituting it along with expressions in (4.6) into (4.5), and separating the real part from the imaginary part, we have which yields where and Notice that for (4.10), it is easy to see that if either P n > 0 and Q n > 0 or P 2 n − 4Q n < 0, then (4.10) has no positive root.Suppose that Q n > 0, P n < 0 and P 2 n − 4Q n > 0, then (4.10) has two positive roots.In addition, if either Q n < 0 or P n < 0 and P 2 n − 4Q n = 0, then (4.10) has only one positive root.
Case 1 It is easy to see that if the conditions (C 0 ) and (C 1 ) : P n > 0 and Q n > 0 or P 2 n − 4Q n < 0, hold, then (4.10) has no positive roots.Hence, by combining with the assumption 2.1, we know that all roots of (4.5) have negative real parts when τ ∈ [0, +∞) under the conditions (C 0 ) and (C 1 ).
In the following, we mainly consider the case of Q n < 0, that is (4.10) has only one positive root ω n .In the following, we will discuss the case which is used to guarantee Q n < 0 under the condition (C 0 ).When τ > 0, according to (4.5) and (4.11), we can define and and then by a simple analysis, we have Γ n (0) = J n (0) > 0 for any n ∈ N 0 .Therefore, the sign of Q n coincides with that of Q n , and in order to guaranteeing Q n < 0, we only need to study the case of Q n < 0. hold, then (4.12) has two positive roots.Without loss of generality, we assume that the two positive roots of (4.12) are x 1 and x 2 , i.e., where Since . By using a geometric argument, we can conclude that where n ∈ N. Therefore, (4.10) has one positive root ω n for n 1 < n < n 2 with n ∈ N, where Furthermore, by combining with the second mathematical expression in (4.9), and notice that a 12 < 0, T n < 0 under the condition (C 0 ), then we have sin(ω n τ ) > 0. Thus, from the first mathematical expression in (4.9), we can set where Re(λ(τ )) represents the real part of λ(τ ).
Proof By differentiating the two sides of with respect to τ , where T n and J n (τ ) are defined by (4.6), we have Therefore, by (4.17), we have Re dλ(τ ) dτ τ =τn,j .
(4.18) Furthermore, according to (4.9), we have Moreover, by combining with (4.18), (4.19) and This, together with the fact that completes the proof, where sign(.)represents the sign function.
Moreover, according to the above analysis, we have the following results.

Numerical simulations
In this section, we verify the analytical results given in the previous sections by some numerical simulations and investigate the direction and stability of Hopf bifurcation.We use the following initial conditions for the system (4.2) u(x, t) = u 0 (x), v(x, t) = v 0 (x), t ∈ [−τ, 0] .For the parameters ℓ = 2, d 11 = 2, d 22 = 3, d 21 = 18, ξ = 0.06, β = 0.5, m = 0.5, s = 0.8, according to Proposition 4.5, we know that system (4.2) undergoes a Hopf bifurcation at τ 1,0 = 6.1498.Furthermore, the direction and stability of Hopf bifurcation can be determined by calculating K 1 K 2 and K 2 using the procedures developed in Section 2. After a direct calculation using MATLAB software, we obtain which implies that the Hopf bifurcation at τ 1,0 = 6.1498 is supercritical and stable.

Conclusion and discussion
In this paper, the diffusive predator-prey system with spatial memory and predator-taxis is proposed, and we derive an algorithm for calculating the normal form of Hopf bifurcation in this system.As a real application, we consider the Holling-Tanner model with spatial memory and predator-taxis.By carrying out the stability and Hopf bifurcation analysis for the Holling-Tanner model with spatial memory and predator-taxis, Yehu Lv the inhomogeneous spatial patterns, i.e., two stable spatially inhomogeneous periodic solutions are found.Furthermore, the supercritical and stable mode-1 and mode-2 Hopf bifurcation periodic solutions are found by using our newly developed algorithm.Numerical simulations verify the results of our theoretical analysis and give us a more intuitive display.
The normal form of Hopf bifurcation for the system (4.2) can be calculated by using our newly developed algorithm in Section 2. Here, we give the detail calculation procedures of B 1 , B 21 , B 22 , B 23 steps by steps.