Turing-Hopf bifurcation of a class of modified Leslie-Gower model with diffusion

In this paper, the dynamics of a class of modified Leslie-Gower model with diffusion is considered. The stability of positive equilibrium and the existence of Turing-Hopf bifurcation are shown by analyzing the distribution of eigenvalues. The normal form on the centre manifold near the Turing-Hopf singularity is derived by using the method of Song et al. Finally, some numerical simulations are carried out to illustrate the analytical results. For spruce budworm model, the dynamics in the neighbourhood of the bifurcation point can be divided into six categories, each of which is exactly demonstrated by the numerical simulations. Then according to this dynamical classification, a stable spatially inhomogeneous periodic solution has been found, which can be used to explain the phenomenon of periodic outbreaks of spruce budworm.

1. Introduction. Since introduced by Leslie and Gower [13,14], the Leslie-Gower (LG) model and its various modifications have received great attention [18,27,7,17,2,3,1,5,23]. Model (1) assumes that the carrying capacity of the prey N is limited by a fixed value A B and the carrying capacity for the predator P is directly proportional to the prey N .
In [19], Murray proposed a class of modified Leslie-Gower model: where R(N ) can be one of the predation terms below Obviously, model (2) can become into the Leslie-Gower model (1) when R(N ) = AN .
In the paper, taking into account the inhomogeneous distribution of the preys and predators in different spatial locations, we consider the following modified Leslie-Gower model with diffusion and Neumann boundary conditions x ∈ Ω, t > 0, where d 1 , d 2 > 0 is the diffusion coefficient characterizing the rate of the spatial dispersion of the prey and predator population, respectively, r, δ > 0 is the linear birth rate of prey and predator, respectively, K > 0 is the carrying capacity of prey population, h > 0 is the proportionality coefficient of prey density to the carrying capacity for the predator, ν is the outward unit normal vector on ∂Ω, Ω = (0, lπ) and the function R satisfies (H 0 ) : R(0) = 0, R(s) is strictly monotone increasing and adequately smooth for s ∈ [0, ∞).
Clearly, the assumption (H 0 ) holds when R(N ) is one of the predation terms in (3). As is well known, Turing-Hopf singularity is a degenerate case where Hopf and Turing bifurcations occur simultaneously. At the moment, the corresponding characteristic equation has a pair of simple purely imaginary roots and a simple zero root. There exist very rich dynamics near Turing-Hopf singularity, which includes stable constant equilibrium, nonconstant steady state, spatially homogeneous and inhomogeneous periodic solutions. Due to the increasing interests in the studying of the bifurcation phenomena in the reaction-diffusion systems arising from the biology, chemistry and physics [4,9,10,11,15,16,21,25,26], and also the fact that the existence of stable spatially inhomogeneous periodic solution can be used to explain the periodic fluctuation of biological populations, we will focus our thoughts on the degenerate case in this paper.
We would also like to mention that Faria presented an approach in reference [6], by which an explicit normal form on the centre manifold near an equilibrium of partial functional differential equations can be calculated, especially for the investigation of generic Hopf bifurcations. However, Faria did not address approaches of calculating normal forms for partial functional differential equations with Turing-Hopf singularity. Although, Song et al. presented a rigorous procedure for calculating the normal form associated with the Turing-Hopf bifurcation of partial functional differential equations in [20]. There are still very few studies on Turing-Hopf bifurcation of the model with practical significance (see [22,24]).
The rest of the paper is organized as follows. In Section 2, we investigate the existence and stability of positive equilibrium and the existence of Turing-Hopf bifurcation. In Section 3, we compute the normal forms on the centre manifold for Turing-Hopf bifurcation by using the method in [20]. In Section 4, we carry out some numerical simulations to support and extend our analytical results.
2. Stability and Turing-Hopf bifurcation. In this section, we consider the stability and Turing-Hopf bifurcation of the following system where f satisfies the assumption (H 0 ). Clearly, system (4) can be written as system (5) by setting The positive equilibrium E * = (u * , v * ) of system (5) satisfies Since the function is strictly monotone decreasing with g(0) = 1 and g(∞) = −∞, we have the following conclusion.
In the following, we analyze the stability and the existence of Turing-Hopf bifurcation for the positive equilibrium E * . Notice that the value of u * only depends on the function f and has nothing to do with the parameters r and δ. Therefore, we take the parameters r and δ as bifurcation parameters. The characteristic equations for the positive equilibrium E * are the following sequence of quadratic equations where N 0 = {0} ∪ N and N is the positive integer set. That is., where with Denote Then, we have that T 0 = δ − α(u * )r and h 0 = β(u * )rδ. Therefore, if α(u * ) ≤ 0, we obtain T 0 > 0 and h 0 > 0 for any δ, r > 0. Which implies that the positive equilibrium E * of system (5) without diffusion is asymptotically stable for any δ, r > 0. If α(u * ) > 0, we know that h 0 > 0 such that for system (5) without diffusion, Hopf bifurcation occurs when T 0 = 0. That is Hopf bifurcation occurs at the line of δ = H 0 (r), where at which Eq.(7) with n = 0 has a pair of purely imaginary roots ±i √ h 0 . It is easy to verify that the following transversality condition holds: Therefore, the result on the system (5) without diffusion follows immediately.
In the following, we focus on the diffusion-driven instability and the bifurcation analysis for system (5) under the case of α(u * ) > 0. For convenience, we define and make some hypotheses as follows: First it is the case for no Turing instability.
Theorem 2.4. If α(u * ) > 0 and either (C 1 ) or (C 2 ) holds, then there is no diffusion-driven Turing instability. In this case, diffusion does not change the stability of the positive equilibrium E * , i.e. the stable and unstable regions are same as Lemma 2.2 (ii).
Proof. It is easy to see from (8) that for any n ∈ N 0 , T n > 0 can be satisfied provided that δ > H 0 (r) and h n = 0 is equivalent to Clearly, if all lines determined by δ = S n (r) lie below the Hopf bifurcation line defined by δ = H 0 (r) for r > 0, then there is no Turing instability. After a simple calculation, we obtain that Obviously, if d 2 ≤ d 1 , we have P n (r) > 0 for any n ∈ N 0 and r > 0, which means that H 0 (r) > S n (r) for any n ∈ N 0 and r > 0. Nextly, we consider the case of d 2 > d 1 and Λ < 0. Firstly, we have S 0 (r) = 0 for any r > 0. This leads to H 0 (r) > S 0 (r) for any r > 0. For n ∈ N, we can obtain that P n (r) > 0 for any r > 0 if Λ < 0. Thus, H 0 (r) > S n (r) for any n ∈ N and r > 0. The proof is complete.
Then it is the case for the existence of Turing instability.
Theorem 2.5. If α(u * ) > 0 and (C 3 ) hold, then there will be diffusion-driven Turing instability and Turing-Hopf bifurcation at some Turing-Hopf singularity. More precisely, denote where p ± n and q + (n, m) are defined below and then we have the following results: (i) The positive equilibrium E * is asymptotically stable for any r > 0, δ > L n (r), n ∈ N 0 . (ii) The critical line of Turing instability is δ = L n (r)(L n (r) > H 0 (r)), n ∈ N and Turing instability occurs for r > 0, H 0 (r) < δ < L n (r), n ∈ N. (iii) System (5) undergoes Turing-Hopf bifurcation at the point (p − 1 , α(u * )p − 1 ); If p − n+1 > p + n , n ∈ N, then system (5) undergoes Turing-Hopf bifurcation at the Proof. From the proof of Theorem 2.4, we know that Turing instability occurs for the condition (C 3 ) holds. In this case, H 0 (r) − S 0 (r) = 0 has no positive root and H 0 (r) − S n (r) = 0, n ∈ N has two positive roots such that It is easy to see that for any m, n ∈ N, Q nm (r) has and only has one positive zero such that S m (r) < S n (r), if 0 < r < q + (n, m), It follows that δ = S m (r) and δ = S n (r) have and only have one intersecting point for any m, n ∈ N, m > n.
From (14) and (15), we know that p ± n increases and tends towards positive infinity with the increase of n and q + (n, m) increases and tends towards positive infinity with the increase of m for a fixed n ∈ N.
Furthermore, we can obtain the transversality condition below According to the above analysis, it is clear that for 0 < r ≤ γ 0 = p − 1 , the boundary of stable region of the positive equilibrium E * is H 0 (r). If p − 2 ≤ p + 1 , then for γ 0 < r ≤ γ 1 = q + (1, 2), the boundary of stable region is S 1 (r) (see Fig.1(a)). Otherwise, when γ 0 < r ≤ γ 1 = p − 2 , the boundaries consist of S 1 (r) with γ 0 < r ≤ p + 1 and H 0 (r) with p + 1 < r ≤ γ 1 (see Fig.1(b)). By using the mathematical induction, we can obtain that the whole boundaries of stable region consist of L n (r), n ∈ N 0 . This, together with the transversality condition (16) completes the proof of (i) and (ii).

XIAOFENG XU AND JUNJIE WEI
Define the real-valued Sobolev space with the inner product Then system (18) can be written as the following differential equation in the abstract where For the formal Taylor expansions of L, It is well known that the eigenvalues of D∆ on X are and the corresponding normalized eigenfunctions are given by where e j is the unit coordinate vector of R 2 and n is often called wave number. Let Then, by a straightforward calculation, we can obtain that ξ 0 ∈ C 2 and ξ n * ∈ R 2 are the eigenvectors associated with the eigenvalues iω c and 0, respectively, and η 0 ∈ C 2 and η n * ∈ R 2 are the corresponding adjoint eigenvectors, where where I 2 is an 2 × 2 identity matrix and Then the phase space X for (19) can be decomposed as where dimP = 3 and π : X → P is the projection defined by 0 ) and β n * = col(β (1) n * , β n * ). According to (23), U ∈ X can be decompose as where (24) can be rewritten as In order to simply the notation, let Then, denoting by L 1 the restriction of L to Q, system (19) is equivalent to a system of abstract ODEs in R 3 × Q, with finite and infinite dimensional variables also separated in the linear term. where and According to [20], by a recursive transformation, we can obtain that the normal form for Turing-Hopf bifurcation reads aṡ where and Clearly, we still need to compute C ijk , D ijk and E ijk . Letting 0, 0, 0, 0) ∂u j1 ∂v j2 , k = 1, 2, j 1 + j 2 = 2, 3.
Then from (20), we have and
4. An example and simulations. In this section, we make some numerical simulations to support and extend our analytical results. Choose R(N ) = AN 2 B+N 2 , then (4) becomes We would like to mention that the system (32) without the diffusion was proposed as a prey-predator model to describe the interaction of spruce budworm and birds, see [19, p71]. In fact, spruce budworm is one of the most destructive insect in North American forests and there are periodic outbreaks of spruce budworm (every 30-40 years lasting for about 10 years) causing billions of dollars loss to forest industry. Therefore, understanding the reason why spruce budworm population can outbreak periodically, is very important to control the growth of spruce budworm and protect spruce and fir forests.
Define the critical bifurcation lines are as follows: H 0 : µ 2 = 0.5152µ 1 ; S 1 : µ 2 = 1.4564µ 1 ; Then, according to the results in [8], the bifurcation diagram in the (µ 1 , µ 2 ) parameter plane and the corresponding phase portraits of system (34) in the (ρ, ς) plane can be shown in Fig.2. The four straight lines H 0 , S 1 , T 1 and T 2 divide the (µ 1 , µ 2 ) parameter plane into six regions marked by x-}. Noticing that the equilibria A 0 , A 1 , A ± 2 and A ± 3 of the normal form system (34) correspond to the constant equilibrium, the spatially homogeneous periodic solution, the nonconstant steady state like cos x-shape and spatially inhomogeneous periodic solution of system (33), respectively. So the dynamics of the original system (33) near the Turing-Hopf point P * in the plane of parameters r and δ can be identified in terms of the dynamics of the normal form system (34). In region x, system (34) has only one equilibrium A 0 and it is asymptotically stable. This leads to the positive constant equilibrium E * of the original system (33) is asymptotically stable, as shown in Fig.3 for (µ 1 , µ 2 ) = (−0.01, 0.02) and the initial value u(x, 0) = 0.1296 + 0.005 cos x, v(x, 0) = 0.0167 + 0.01 cos x.
In region y, system (34) has three equilibria: A 0 , A + 2 and A − 2 . The equilibrium A 0 is unstable and the equilibria A + 2 and A − 2 are asymptotically stable. So, the original system (33) has an unstable constant equilibrium and two stable nonconstant steady states like cos x-shape. For the fixed parameter (µ 1 , µ 2 ) = (0.022, 0.014) and choosing different initial values, system (33) can converge to one of these two nonconstant steady states, as shown in Fig.4(a)     In region z, system (34) has four equilibria: A 0 , A 1 , A + 2 and A − 2 . The equilibria A 0 and A 1 are unstable and the equilibria A + 2 and A − 2 are still asymptotically stable. So, the original system (33) has an unstable constant equilibrium, an unstable spatially homogeneous periodic solution and two stable nonconstant steady states like cos x-shape. Choosing (µ 1 , µ 2 ) = (0.02, 0.01) and the initial value as u(x, 0) = 0.1576 − 0.002 cos x, v(x, 0) = 0.0234, the dynamics of system (33) evolves from the spatially homogeneous periodic solution, and spatially inhomogeneous periodic solution and finally to the nonconstant steady state, as shown in Fig.5.