Double Hopf bifurcation of a diffusive predator–prey system with strong Allee effect and two delays

Abstract. In this paper, we consider a diffusive predator–prey system with strong Allee effect and two delays. First, we explore the stability region of the positive constant steady state by calculating the stability switching curves. Then we derive the Hopf and double Hopf bifurcation theorem via the crossing directions of the stability switching curves. Moreover, we calculate the normal forms near the double Hopf singularities by taking two delays as parameters. We carry out some numerical simulations for illustrating the theoretical results. Both theoretical analysis and numerical simulation show that the system near double Hopf singularity has rich dynamics, including stable spatially homogeneous and inhomogeneous periodic solutions. Finally, we evaluate the influence of two parameters on the existence of double Hopf bifurcation.


Introduction
In nature, population growth is limited by environmental resources, this explains why the population cannot grow indefinitely. A proper population growth rate is crucial in describing the population dynamics. The logistic growth proposed by Pearl and Reed in 1920 [18] has been studied by many scholars since then. It assumes that per capita growth rate is a monotonically decreasing function of the population density. The logistic growth considers the intraspecific competition among individuals for limited resources. In ecology, intraspecific competition involves competition for food, mates, shelters, and so on. However, for many social animals, the individual can benefit from the presence of conspecifics, the intraspecific cooperation cannot be ignored. Behaviours of intraspecific cooperation include cooperative predation, cooperative defense against enemies, etc.
The strong Allee effect growth characterizes the population growth of both intraspecific cooperation and intraspecific competition [1,8,13], which means the growth rate would be negative when the population density is below a threshold [15,19]. For the strong Allee effect, the most commonly used expression is where 0 < β < K is the Allee threshold. Obviously, when u < β, the growth is negative, and the population will eventually die out. This phenomenon is abundant in ecology [4,13], include inbreeding depression, failure to satiate natural enemies, failure in mate finding, and so on.
In population dynamics, delay occurs in almost every situation, and its effect cannot be neglected [14]. Time delays tend to destabilize the system and lead to oscillatory behavior. Indeed, population densities of many species are known to fluctuate nearly periodically over time [20], a phenomenon to which the delay may provide an explanation. Notice that both the positive feedback of intraspecific cooperation and the negative feedback of intraspecific competition can have delay to the growth of the population. Chang et al. [5] investigated the dynamics of a scalar population model with delayed Allee effect as follows: They obtained rigorous stability and Hopf bifurcation results for the above model. They showed that the increasing of τ may enlarge the basin of attraction of u = 0 and lessen the basin of attraction of u = 1. Besides, they showed that large delay may devote to the extinction of the population. For many species, the delay feedback in intraspecific cooperation is not necessarily the same as that in intraspecific competition. Jankovic and Petrovskii investigated two models as follows [13]:U They found that the inclusion of delay in intraspecific cooperation term does not necessarily confer instability, but a delay in intraspecific competition always results in instability, thus leads to population cycles or even extinction. Besides, they obtained that delay in competition term has more dominant effects on population dynamics compared to the delay in cooperation term.
Besides, the digestion or gestation delay plays an important role in characterizing the increase of the predator growth. Recently, many researchers consider the digestion delay in their models [3,6,7,17,23]. Chen et al. [6] proposed a diffusive predator-prey model with digestion delay. They investigated the effect of the digestion delay on the system and obtained the stability of the equilibria and the existence of Hopf bifurcation.
Motivated by the previous works, we introduce a delay in intraspecific competition in prey, consider the digestion delay in predator, and conclude a diffusive predator-prey system with strong Allee effect with two delays as follows: where u = u(x, t), v = v(x, t) denote the density of the prey and predator at location x and time t, respectively; d 1 , d 2 mean the diffusion coefficients; τ 1 , τ 2 denote the delay feedback of intraspecific competition in prey and the digestion delay in predator, respectively; c denotes the Allee threshold, and K is the carrying capacity; m represents the modified capture rate; b denotes the semi-saturation coefficient, and d the death rate of predator. The homogeneous Neumann boundary condition is imposed so that there is no population flow across the boundary, ν denotes the outward normal to the boundary ∂Ω. For simplicity, we choose Ω as a spatial interval [0, lπ].
Our purpose is to explore the joint effect of two delays on the dynamics of system (1). The introduction of one delay may unstabilize the equilibria, lead to the occurrence of Hopf bifurcation [24]. The dynamics of the system with two delays would be more complex. There are several related studies on systems with two delays [2,21,22]. Most of the works regard one delay as a constant and the other as variable. To evaluate the joint effect of two delays, it is necessary to take two delays as variables. Gu et al. [11] and Wang et al. [16] provide methods in studying the stability of a system with two delays vary simultaneously. Du et al. investigated a diffusive Leslie-Gower model with two delays [10], they performed the stability analysis and explored the existence of double Hopf bifurcation.
The highlights of this paper are as follows. For one thing, the Hopf bifurcation curves and double Hopf bifurcation curves on a two-parameter plane are given. For another, the calculation of the norm forms near the double Hopf singularity need subtle analysis, since two simultaneously varying delays make the calculation process more complex. Besides, the stable spatially homogeneous and inhomogeneous periodic solutions are theoretically proved and illustrated near the double Hopf singularities.
Our paper is organized as follows. In Section 2, we obtain the stability switching curves of the positive equilibrium and present the Hopf and double Hopf bifurcation theorem. In Section 3, we calculate the normal forms on the center manifold near the double Hopf singularity. In Section 4, we carry out some numerical simulations for illustrating the theoretical results. Finally, conclusions are made to sum up the paper.

Stability switching curves and the existence of double Hopf bifurcation
In this section, we perform the stability and bifurcation analysis of the system. http://www.journals.vu.lt/nonlinear-analysis

Stability switching curves
By direct calculation, system (1) has three constant steady states E 0 (0, 0), E 1 (c, 0), E 2 (K, 0) and a possible positive constant steady state To guarantee the existence of E * , we make the following assumption: Since the study of boundary equilibria is of little biological significance, in this paper, we focus on the study of E * and assume that (H0) is always true. The linearization of system (1) at E * is Then the characteristic equation of E * is where I 2 is the 2 × 2 identity matrix, and M n = −(n 2 /l 2 )D. For simplicity, denote Then (2) is equivalent to D n (λ; τ 1 , τ 2 ) := P 0,n (λ) + P 1,n (λ)e −λτ1 + P 2,n (λ)e −λτ2 = 0, n ∈ N 0 , ( where To maintain the stability of E * at τ 1 = τ 2 = 0, we further make the hypothesis Then we can easily get the following theorem. Theorem 1. If (H0) and (H1) hold, E * is locally asymptotically stable for τ 1 = τ 2 = 0. Now we will investigate the joint effect of τ 1 and τ 2 on the stability of E * under (H0) and (H1). To apply the method of the stability switching curves, we first check assumptions (I)-(IV) in [16] for all n ∈ N 0 . Obviously, (I)-(IV) are true for P 0,n (λ), P 1,n (λ) and P 2,n (λ) in Eq. (4). Particularly, assumption (IV) maintains feasible ω s to be bounded.
To proceed the analysis, we introduce the following lemma [21].
Lemma 1. As (τ 1 , τ 2 ) varies continuously in R 2 + , the number of zeros (with multiplicity counted) of D n (λ; τ 1 , τ 2 ) on C + can change only if a zero appears or cross the imaginary axis.
Lemma 2. Ω n consists of a finite number of intervals of finite length.
The above lemma can be easily proved, here we omit the proof. In the following, we will investigate the relationship between different stability switching curves. From Eqs. (8) and (9), for any ω ∈ Ω n , on the corresponding stability switching curves, one can obtain either τ + 2,j2,n or τ − 2,j2,n for a given τ + 1,j1,n . With the aid of MATLAB, we verify that when τ 1 = τ + 1,j1,n , τ 2 = τ − 2,j2,n ; similarly, when τ 1 = τ − 1,j1,n , τ 2 = τ + 2,j2,n . Therefore, Proposition 1. T n is the set of all stability crossing curves on τ 1 ,τ 2 -plane for D n (λ; . By direct calculation, we can get that Obviously, (11) means T +k j1,j2,n connects with T −k j1+δ a 1 ,j2−δ a 2 ,n and T −k 2 ,n at its two ends. In (11) , one should in mind that the values of δ a k,n i or δ b k,n i might be different for different k, n and i.

Crossing directions
To present the Hopf bifurcation and double Hopf bifurcation theorem, we define the crossing directions of the points on the stability switching curves in this subsection.
Let λ = σ +iω. By the implicit function theorem, τ 1 , τ 2 can be expressed as functions of σ and ω under some non-singular condition. For convenience, denote and for l = 1 and 2, We can easily verify that Since T ±k j1,j2,n are piecewise differentiable, by the implicit function theory, we have as long as the nonsingular condition R 1 I 2 − R 2 I 1 = 0 holds. For any crossing curve T ±k j1,j2,n , we define the direction of the curve corresponding to the increasing ω ∈ Ω n,k the positive direction. The region on the left-hand (righthand) side when we move in the positive direction of the curve the region on the left (right). By direct calculation, we get that the tangent vector of T ±k j1,j2,n along the positive direction is (∂τ 1 /∂ω, ∂τ 2 /∂ω), the normal vector of T ±k j1,j2,n pointing to the right region is (∂τ 2 /∂ω, −∂τ 1 /∂ω). Besides, a pair of complex characteristic roots across the imaginary axis from the left to the right of the complex plane as σ increases from negative to positive, thus (τ 1 , τ 2 ) moves along (∂τ 1 /∂σ, ∂τ 2 /∂σ). We can deduce that if then the region on the right of T ±k j1,j2,n has two more characteristic roots with positive real parts. If δ(ω) < 0, then the left region of T ±k j1,j2,n has two more characteristic roots with positive real parts. Since det R0 −I0 One can also verify that Hence, sign δ(ω ∈Ω k,n ) = ± sign(ω 2 |P 2P3 − P 0P1 | sin(ψ 1,n )), whereΩ k,n denotes the interior of Ω k,n .
Lemma 3. For any k = 1, 2 . . . , r n , we have Therefore, the region on the right (left) of the crossing curve T +k j1,j2,n (T −k j1,j2,n ) has two more characteristic roots with positive real parts.

Hopf and double Hopf bifurcation theorem
Now we introduce the Hopf bifurcation theorem [10] and conclude the double Hopf bifurcation theorem.
Theorem 2. For any j = 1, 2 . . . , r n , T j n is a Hopf bifurcation curve in the following sense: for any p ∈ T j n and for any smooth curve Γ intersecting with T j n transversely at p, we define the tangent of Γ at p by l. If ∂ Re λ/∂ l| p = 0 and the other eigenvalues of (3) at p have nonzero real parts, then system (1) undergoes a Hopf bifurcation at p when parameters (τ 1 , τ 2 ) cross T j n at p along Γ .
If there exist ω j1,k1 and ω j2,k2 such that T j1 k1 and T j2 k2 intersect, then there are two pairs of pure imaginary roots of (3) at the intersection. Thus, system (1) may undergoes double Hopf bifurcation at this intersection under certain conditions. Here we summarize the following double Hopf bifurcation theorem.
Theorem 3. If two Hopf bifurcation curves T j1 k1 and T j2 k2 intersect at q, we define the tangent of T j1 k1 and T j2 k2 at q by l 1 and l 2 , respectively. In addition, if l 1 and l 2 are linearly independant, namely, for α 1 , α 2 ∈ R, α 1 l 1 + α 2 l 2 = 0 can only arrive at α 1 = α 2 = 0, then system (1) undergoes a double Hopf bifurcation at q. Now we can summarize the following theorem about the stability of E * .
This theorem can be easily proved with the aid of Lemma 1. For a better understanding of this theorem, please refer to Fig. 7(a).

Normal form of double Hopf bifurcation
From Section 2 we get that system (1) may undergo double Hopf bifurcations near E * (u * , v * ). To investigate the dynamics near the double Hopf singularity, we will calculate the normal forms on the center manifold near the double Hopf singularity by the method derived in [9].
When P 5 (0.84, 7.95) falls in region R 1 , E * is a sink. When P 6 (0.84, 7.9662) is chosen in R 2 or P 7 (0.84, 7.9692) in R 3 , there exists a stable periodic solution. We only draw the case for P 7 ∈ R 3 in Fig. 10.     Since the strong Allee effect is considered in system (1), the over-predation phenomenon will occur when the initial value of the predator is large enough compared to the prey. We only take one example to show the over-predation phenomenon, which is revealed in Fig. 11. Now we use parameters in (21) as a baseline to show the sensitivity of c and K on model (1) by numerical simulations.
First we consider the sensitivity of c on the stable region of E * . We draw the left-most and bottom stability crossing curves and zoom them in [0.65, 1] × [0, 14] in Figs. 12(a) and (b). Then K is chosen as the perturbation parameter with all other parameters fixed as the same in (21). We also draw the left-most and bottom stability crossing curves and zoom them in [0.65, 1] × [0, 16] in Figs. 12(c) and (d).
From Fig. 12 we can conclude that although the influence of a slight perturbation of c or K on the stability region cannot be evaluated quantitatively, the double Hopf bifurcation still exist, which indicates that the double Hopf bifurcation in the system has strong anti-interference ability.

Conclusion
A diffusive predator-prey system with two delays and strong Allee effect is investigated in our paper. This system considers both the delay feedback of intraspecific competition in prey and the delay feedback of digestion in the predator. By applying the method of stability switching curves, we study the joint effect of the two delays on the stability of the positive constant steady state. With the aid of the crossing directions, we present the Hopf bifurcation theorem and give sufficient conditions for the existence of the double Hopf bifurcation.
To explore the dynamics of the system near the double Hopf singularities, we calculate the normal forms on the center manifold near the double Hopf singularities. The calculation formulas of the normal form we give here are at a double Hopf singularity with k 1 = 0, k 2 = 1. Then we get the corresponding unfoldings and the bifurcation sets. With a set of parameters, we theoretically prove and illustrate the existence of homogeneous and inhomogeneous periodic solution. Besides, heteroclinic orbits are found near the double Hopf singularity.
Populations with strong Allee effect can be wiped out by the over-predation phenomenon, we prove this phenomenon numerically. Finally, we carry out the sensitivity analysis to show the impact of two parameters on the dynamics of the system. Numerical simulation shows that small changes of c or K normally won't affect the existence of double Hopf singularities of the system.