1. 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 intraspe- cific 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 . = 0 and lessen the basin of attraction of . = 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]:

They found that the inclusion of delay in intraspecific cooperation term does not necessar- ily 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 . = .(x, t), . = .(x, t) denote the density of the prey and predator at location . and time ., respectively; .., .. mean the diffusion coefficients; .., .. denote the delay feedback of intraspecific competition in prey and the digestion delay in predator, respectively; . denotes the Allee threshold, and . is the carrying capacity; . represents the modified capture rate; . denotes the semi-saturation coefficient, and . 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.

2. Stability switching curves and the existence of double Hopf bifurcation

In this section, we perform the stability and bifurcation analysis of the system.

2.1 Stability switching curves

By direct calculation, system (1) has three constant steady states ..(0, 0), ..(c, 0), ..(K, 0) and a possible positive constant steady state ..(.., v.), where .. = bd/(.−.), .. = .(.. .)(.. + .)(1 ../K./m. To guarantee the existence of .., 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 .. and assume that (H0) is always true.

The linearization of system (1) at .. is

where

Then the characteristic equation of .. is

where .. is the 2 × 2 identity matrix, and M. = −(../l.).. For simplicity, denote

Then (2) is equivalent to

where

To maintain the stability of .. at .. = .. = 0, we further make the hypothesis (H1) (. − ..)(2.. + . − .).(. + ..) < u. − ..

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 N0. Obviously, (I)–(IV) are true for P0,n(λ), P1,n(λ) and P2,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 (.., τ.) varies continuously in R2 , the number of zeros .with multiplicity counted. of D.(.; .., τ.) on C. can change only if a zero appears or cross the imaginary axis.

Since . = 0 is not a zero of (3), we will seek purely imaginary zeros of (3) to explore the stability switching curves. Assume . = i. (ω > 0) is a zero of (3), then

ince |e^−iωτ2 | = 1, after simplification, we obtain

where ._1,n(.) = Re(−._0,n(i.). _1,n(i.)), ._1,n(.) = Im(−._0,n(i.). _1,n(i.)).

Without loss of generality, assume A21;n(!) + B21;n(!) > 0, then there exists somecontinuous function -1;n(!) 2 Argf-P0;n(i!)P1;n(i!)g \ (-π; 2π] such that (5) becomes

then the existence of .. ∈ R. in (6) equals

Obviously, Eq. (7) contains the special case ..(.) + ..(.) = 0. We denote .. to be

the collection of positive .’s such that (7) hold. Let ._1,n(.) ∈ [0, π] satisfies

then we have

Similarly, we can get that

where . satisfies the following condition:

Denote .. as the collection of all .’s such that (10) hold, by squaring both side of (7) and (10), we can get that .. = ... For simplicity, denote Ω. := .. = ...

To clarify that the stability switching curves of a finite region contains all the stability switching curves in that region, we present the following lemma.

Lemma 2

Ω. 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 switch- ing curves. From Eqs. (8) and (9), for any . ∈ Ω., on the corresponding stability

And

Proposition 1

T. is the set of all stability crossing curves on τ.,τ.-plane for D.(.; .., τ.) = 0. For any (.., τ.) ∈ T., D.(.; .., τ.) = 0 has at least one root .ω with ω . Ω..

Notice that F..a_k,n) = F..b_k,n) = 0, then ψ_i,n.a_k,n) = δ.k,n π, ψ_i,n.b_k,n) =

By direct calculation, we can get that

11Obviously, (11) means T +kj1;j2;n connects with T 􀀀kj1+a1 ;j2􀀀a2 ;n and T 􀀀kj1+b1;j2􀀀b2;n at its twoends. In (11), we denote ai ; bi instead of ak;ni ; bk;ni , one should in mind that the valuesof ak;ni or bk;ni might be different for different k, n and i.

2.2 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, .., .. 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

are piecewise differentiable, by the implicit function theory, we haveParrafo con ecuacion

as long as the nonsingular condition R1I2 􀀀 R2I1 6= 0 holds.For any crossing curve T kj1;j2;n, we define the direction of the curve correspondingto the increasing ! 2 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 kj1;j2;n along the positivedirection is (@1=@!; @2=@!), the normal vector of T kj1;j2;n pointing to the right region is(@2=@!;􀀀@1=@!). Besides, a pair of complex characteristic roots across the imaginaryaxis 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 has two more characteristic roots with positive realparts. If (!) < 0, then the left region of has two more characteristic roots with

positive real parts. Since det􀀀 R0 􀀀I0I0 R0= R20 + I20 > 0, we get that if R20 + I20 6= 0,sign (!) = signfR1I2 􀀀 R2I1g. One can also verify that

Hence, sign .(ω Ω.k,n) = sign(.. ...¯3 ...¯1 sin(.._,n)), where .˚k,n denotes the interior of Ω_k,n.

Lemma 3

For any k = 1, 2 . . . , r., we have

PARRAFO CON ECUACIONTherefore, the region on the right (left) of the crossing curve () has twomore characteristic roots with positive real parts.

2.3 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 . . . , rn, T j is a Hopf bifurcation curve in the following sense: for any p ∈ T j and for any smooth curve Γ intersecting with T j transversely at p,For any j = 1; 2 : : : ; rn, T jn is a Hopf bifurcation curve in the followingsense: for any p 2 T jn and for any smooth curve 􀀀 intersecting with T jn transversely at p,we define the tangent of 􀀀 at p by~l. If @ Re =@~ljp 6= 0 and the other eigenvalues of (3)at p have nonzero real parts, then system (1) undergoes a Hopf bifurcation at p whenparameters (1; 2) cross T jn at p along 􀀀.If there exist !j1;k1 and !j2;k2 such that T j1k1 and T j2k2 intersect, then there are twopairs of pure imaginary roots of (3) at the intersection. Thus, system (1) may undergoesdouble Hopf bifurcation at this intersection under certain conditions. Here we summarizethe following double Hopf bifurcation theorem.

Theorem 3

If two Hopf bifurcation curves T j1k1and T j2k2intersect at q, we define thetangent of T j1k1and T j2k2at q by~l1 and~l2, respectively. In addition, if~l1 and~l2 are linearlyindependant, namely, for 1; 2 2 R, 1~l1 + 2~l2 = 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.

Theorem 4

Assume that (H0) and (H1) are true. For any point P˜(τ¯1, τ¯2) on the τ1,τ2- plane, if there exists a curve segment ˜l connecting P˜ and the origin such that ˜l does not intersect any stability switching curves, then E∗ is stable when τ1 = τ¯1, τ2 = τ¯2.

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).

3. Normal form of double Hopf bifurcation

From Section 2 we get that system (1) may undergo double Hopf bifurcations near ..(.., v.). To investigate the dynamics near the double Hopf singularity, we will cal- culate the normal forms on the center manifold near the double Hopf singularity by the method derived in [9].

Assume that .. > τ., let .¯(x, t) = .(x, τ..) .., .¯(x, t) = .(x, τ..) .. and drop the bars. System (1) can be written as

Where

And

Define the real-valued Hilbert space

and the complexfication space of . by .. := . ⊕ i. = {.. + i.., U., U. ∈ .}

denote the phase spacewith the sup norm. We write

Denote the double Hopf singularity by ( 1 ; 2 ), introduce two parameters = (1; 2)with 1 = 1 􀀀 1 and 2 = 2 􀀀 2 , denote U(t) = [u(t); v(t)]T, then system (12) canbe written as

where

Linearizing system (13) at (0, 0), we have

It is well known that the eigenvalues of .∆ on . are −..../l. and −..../l. with cor

responding normalized eigenfunctions .⁽¹⁾(.) = γ.(.)(1, 0)., .⁽²⁾(.) = γ.(.)(0, 1). with γ.(.) = cos((n/l).).ǁ cos((n/l).)ǁ.2 . Define the subspace of C as B., where B. System (13) can be written as

Where

Formulating (15) in the extended Banach space:

rewrite (15) as an abstract ordinary differential equations on BC, then

where is the infinitesimal generator of the ..-semigroup of solution maps of the linear equation (14) defined by

with dom( ) = . : .˙ , φ(0) dom(∆) with ..(.) = 0 for . [ 1, 0) and ..(0) = .. On ., system (14) is equivalent to the retarded functional differential equation

R2. Thus, there exist matrix functions m 2 BV([􀀀1; 0];R) such that 􀀀(k2m=l2) D0'(0) + L0(') =R 0􀀀1 dm()'(). Let An denote the infinitesimal generator of thesemigroup generated by (17), An denote the formal adjoint of An under the bilinear form

System (14) has two pairs of pure imaginary eigenvalues i...2∗ and i...2∗ at the double Hopf singularity and the other eigenvalues with nonzero real parts. Here we do not consider the strongly resonant cases, i.e., .. : .. =/ . : . for p, q ∈ N and 1 ≤ p, q ≤ 3. Let ..(.) = (1, p₁₂).e^i.1.2∗ . , ..(.) = (1, p₃₂). e^i.2.2∗ . (. ∈ [−1, 0]) are the eigenvectors of A. (. = 1, 2) corresponding to the eigenvalue i...2∗, i...2∗, respectively. Let ..(.) = ..(1, p.12)e⁻ⁱ.1.2∗ s, ψ.(.) = ..(1, p.32)e⁻ⁱ.2.2∗ . (. ∈ [0, 1]) are the eigenvectors of A∗. (. = 1, 2) corresponding to the eigenvalue i...2∗, i...2∗, respectively. Choose P. and ... as the basis of the generalized eigenspace of A., A∗.

corresponding the the eigenvalues {±i...2∗, ±i...2∗}, respectively, which

with and satisfy

A.Φ. = Φ.B.1 , A.. Ψ. = B.. Ψ., (Ψ., Φ... = . with B.1 = diag(i...2∗, −i...2∗)

and B.2 = diag(i...2∗, −i...2∗). By direct calculations, we get that

Now decompose BC into a direct sum of center subspace and its compleΣmentary space:

is the projection defined by

where A1 is the restriction of A on Q1 _ Ker _ ! Ker _; A1_ = A_ for _ 2 Q1.

Consider the formal Taylor expansion .(σ, ϕ) = ..=2(1/j!)G..σ, ϕ), where G. is the .th Fréchet derivation of .. Then we can rewrite (18) as

where . = (.., z., z., z.). ∈ C4, . ∈ Q1, and f. = (. ., f .), . ≥ 2, are given by

For a normed space Y , V ⁴⁺²(. ) denotes the space of homogeneous polynomials of degree . in . ∈ C4, . = (.., σ.) with coefficients in . , we have

and the norm

By a recursive transformations of variables (z, w, σ) = (.ˆ, w., σ) + (1/j!)(. .(.ˆ, σ., U .(.ˆ, σ), 0) with U. = (. ., U .) ∈ . ⁴⁺²(C4) × . ⁴⁺²(Q1), we can obtain the normal forms of (19). We conclude that this recursive process transforms (19) into the following equation:

where g. = (.., g.) have the following form g..z, w, σ) = fj .z, w, σ) − M.U..z, σ),

and Uj 2 V 4+2j (C4) V 4+2j (Q1) satisfying Uj(z; ) = (Mj)􀀀1 ProjIm(M1j )Im(M2j ) fj(z; 0; ), where fj = (f1j ; f2j ) stand for the terms of order j in (z;w), which areobtained after the computation of normal forms up to order j 􀀀 1. The normal formtruncated to the third order has the form

Here The calculationsof g12(z; 0; _) and g13(z; 0; 0) are the same as the calculation in [10]. After the calculations, the normal form truncated to the third order on the center manifold for double Hopf bifurcation can be obtained. By making the polar coordinate transformation calculations, the normal form truncated to the third order on the center manifold for double Hopf bifurcation can be obtained. By making the polar coordinate transformation

where .., ρ. > 0. Take .. = ..√|.₂₁₀₀|, .. = ..√|.₀₀₂₁|, .. = sign(Re .₂₁₀₀), and drop the hats. We obtain the simplified system

where .. = ..(Re .₁₁.. + Re .₂₁..), .. = ..(Re .₁₃.. + Re .₂₃..), .. = .... Re .₁₀₁₁. Re .₀₀₂₁, .. = Re .₁₁₁₀.Re .₂₁₀₀ and .. = ..... As was discussed in [12, Chap. 7.5], there are 12 distinct kinds of unfoldings for (20). In Section 4.1, three cases arise, thus we will draw the bifurcation sets and phase portraits for the unfoldings of case VIII and VIb for illustration.

4. Numerical simulations

In this section, we will do some numerical simulations to check the previous theoretical results. Besides, we will qualitatively analyze the influence of the small perturbations of . and . on the existence of the double Hopf bifurcations.

Now we carry out some simulations for system (1). Choose

One can verify that (H0) holds and .. is (2.5, 0.75). Since .₁₁ + .₁₁ = −0.8333 < 0, (H1) holds, and .. is locally stable for .. = .. = 0.

When .., .. vary simultaneously, we follow the process presented in Section 2.1. As is shown in Fig. 1(a), ..(0) > 0 and ..(.) = 0 has four zeros ._1,0 = 0.06146, ._1,0 = 0.16390, ._2,0 = 1.22056, ._2,0 = 1.44592, the crossing set is ._{1,0 .2,0} =

[._1,0, b_1,0] [._2,0, b_2,0]. For the two ends of ._1,0, we have ._1,0(._1,0) = π, ψ_2,0(._1,0) = 0, ._1,0(._1,0) = 0, and ._2,0(._1,0) = ., thus δ. = 1, δ. = 0, δ. = 0, δ. = 1. From

Eq. (11) we get that T +0;j1;j2 connect T 􀀀0;j1+1;j2 at the left point a1;0 and T 􀀀0;j1;j2􀀀1 at theright point b1;0. As to 2;0, we get 1;0(a2;0) = 0, 2;0(a2;0) = , 1;0(b2;0) = 0, and2;0(b2;0) = 0, thus a1 = 0, b1= 1, a2 = 0, b2= 0. From (11) we get that T +0;j1;j2connect T 􀀀0;j1;j2􀀀1 at a2;0 and T 􀀀0;j1;j2 at b2;0. From Lemma 2 we obtain all the stabilityswitching curves of n = 0 given by T 10 and T 20 , which are showed in Figs. 1(b) and (c),respectively. To give a clearer description of the stability crossing curves, we take the3.5

lower left portion of T 1 and the left-most curve of T 2 as examples and draw the detailed figures in Fig. 2. The numerical results in Fig. 2 coincide with the theoretical analysis.

When . = 1, ..(.) = 0 has three zeros ._1,1 = 0.08886, ._2,1 = 1.28106, ._2,1 = 1.48906, and the crossing set is ._1,1 ∪ ._2,1 = [0, b_1,1] ∪ [._2,1, b_2,1]. Then we get the stability switching curves T 1 and T 2, which are shown in Figs. 3(b) and (c).

When . = 2, ..(.) = 0 has two zeros ._1,2 = 1.39120, ._1,2 = 1.56591, thus ._1,2 = [._1,2, b_1,2]. The stability switching curves . are shown in Fig. 5(a). When . = 3, ..(.) = 0 has two zeros ._1,3 = 1.3822, ._1,3 = 1.5362, thus ._1,3 = [._1,3, b_1,3], the stability switching curves . are shown in Fig. 5(b). When . = 4, ..(.) = 0 has two zeros ._1,4 = 0.9052, ._1,4 = 1.1001, thus ._1,4 = [._1,4, b_1,4], the stability switching curves . are shown in Fig. 5(c). When . ≥ 5, F.(.) > 0 for all . ≥ 0, thus there are no stability switching curves on the ..,τ.-plane.

Combining the stability switching curves shown in Figs. 1, 3, and 5 together, zooming them in one figure in [0, 1.3]×[0, 14], we get the stability switching curves shown in Fig. 6. By Lemma 2, we have plotted all the stability switching curves in [0, 1.3] × [0, 14].

We focus on the lower bottom region bounded by left-most curve of T 1, T 1 and the

lowest curve of . in Fig. 6. By Lemma 1, .. is stable in the bottom left region, since the crossing directions of the three stability switching curves are all pointing outside of this region. These three stability crossing curves intersect at HH.(0.7717,9.8717), HH.(0.8413,7.9598), and HH.(0.9011,7.6982). By the normal form derivation process in Section 3, we obtain the parameters in Table 1. For HH., we get that case VIII oc- curs [12]. Near the double Hopf singularity HH., there are six different kinds of phase diagrams in six regions, which are shown in Fig. 7(a). The six regions are divided by .., . . . , l. with ..: .. = 1.5379.., ..: .. = 62.6720.., ..: .. = 58.5472.. (.. ≥ 0), ..: .. = 15.1454.. (.. ≥ 0). For HH., we get that case VIb occurs. There are also six different kinds of phase diagrams in six regions near the double Hopf bifurcation

point HH., which are shown in Fig. 7(b). The six regions are divided by .., . . . , k. with ..: .. = 3.7480.., ..: .. = 8.0660.., ..: .. = 7.5157.. (.. ≤ 0), ..: .. = 7.3311.. (.. ≤ 0).

When (.., τ.) are chosen as ..(0.72, 9.7) in region .., the positive equilibrium is a sink, which is shown in Fig. 8. When (.., τ.) are chosen as ..(0.775,9.7) in .. or ..(0.775,9.87) in .., or ..(0.775,9.7) in .., there exists a stable periodic solution originating from the Hopf bifurcation. We only draw the case for .. ∈ .. in Fig. 9.

When ..(0.84,7.95) falls in region .., .. is a sink. When ..(0.84,7.9662) is chosen in .. or ..(0.84,7.9692) in .., there exists a stable periodic solution. We only draw the case for .. ∈ .. in Fig. 10.

Since the strong Allee effect is considered in system (1), the over-predation phe- nomenon 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 . and . on model (1) by numerical simulations.

First we consider the sensitivity of . on the stable region of ... 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 . 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).

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.

5 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 cal- culate 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 k1 = 0, k2 = 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 phe- nomenon, 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.

Acknowledgments

The authors thank the two anonymous referees for their very helpful comments which greatly improve the manuscript.

References

1. W.C. Allee, E.S. Bowen, Studies in animal aggregations: Mass protection against colloidal silver among goldfishes, J. Exp. Zool., 61(2):185–207, 1932, https://doi.org/10. 1002/jez.1400610202.

2. Q. An, E. Beretta, Y. Kuang, C. Wang, H. Wang, Geometric stability switch criteria in delay differential equations with two delays and delay dependent parameters, J. Differ. Equations, 266(11):7073–7100, 2019, https://doi.org/10.1016/j.jde.2018.11.025.

3. E. Beretta, Y. Kuang, Convergence results in a well-known delayed predator–prey system, J. Math. Anal. Appl., 204(3):840–853, 1996, https://doi.org/10.1006/jmaa. 1996.0471.

4. D.S. Boukal, L. Berec, Single-species models of the Allee effect: extinction boundaries, sex ratios and mate encounters, J. Theor. Biol., 218(3):375–394, 2002, https://doi.org/ 10.1006/jtbi.2002.3084.

5. X. Chang, J. Shi, J. Zhang, Dynamics of a scalar population model with delayed Allee effect, Int. J. Bifurcation Chaos, 28(12):1850153, 2018, https://doi.org/10.1142/ S0218127418501535.

6. S. Chen, J. Shi, J. Wei, The effect of delay on a diffusive predator–prey system with Holling type-II predator functional response, Commun. Pure Appl. Anal., 12(1):481–501, 2013, https://doi.org/10.3934/cpaa.2013.12.481.

7. K.L. Cooke, Z. Grossman, Discrete delays, distributed delay and stability switches, J. Math. Anal. Appl., 86(2):592–627, 1982, https://doi.org/10.1016/0022-247X(82) 90243-8.

8. F. Courchamp, L. Berec, J. Gascoigne, Allee Effect in Ecology and Conservation, Oxford Univ. Press, 2008, https://doi.org/10.1093/acprof:oso/9780198570301. 001.0001.

9. Y. Du, B. Niu, Y. Guo, J. Wei, Double Hopf bifurcation in delayed reaction–diffusion system, J. Dyn. Differ. Equations, 32(1):313–358, 2020, https://doi.org/10.1007/ s10884-018-9725-4.

10. Y. Du, B. Niu, J. Wei, Two delays induce Hopf bifurcation and double Hopf bifurcation in a diffusive Leslie–Gower predator–prey system, Chaos, 29(1):013101, 2019, https:// doi.org/10.1063/1.5078814.

11. K. Gu, S. Niculescu, J. Chen, On stability crossing curves for general systems with two delays, J. Math. Anal. Appl., 311:231–253, 2005, https://doi.org/10.1016/j. jmaa.2005.02.034.

12. J. Guckenheimer, P. Holmes, Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields, Springer, New York, 1983.

13. M. Jankovic, S. Petrovskii, Are time delay always destabilizing? Revisiting the role of time delays and the Allee effect, Theor. Ecol., .(4):335–349, 2014, https://doi.org/10. 1007/s12080-014-0222-z.

14. Y. Kuang, Delay Differential Equations with Applications in Population Dynamics, Academic Press, New York, 1993.

15. M.A. Lewis, S.V. Petrovskii, J.R. Potts, The Mathematics Behind Biological Invasions, Springer, Cham, 2016, https://doi.org/10.1007/978-3-319-32043-4.

16. X. Lin, H. Wang, Stability analysis of delay differential equations with two discrete delays, Can. Appl. Math. Q., 20(4):519–533, 2012.

17. J. Liu, X. Zhang, Stability and Hopf bifurcation of a delayed reaction–diffusion predator–prey model with anti-predator behaviour, Nonlinear Anal. Model. Control, 24(3):387–406, 2019, https://doi.org/10.15388/NA.2019.3.5.

18. R. Pearl, L. Reed, On the rate of growth of the population of the United States since 1790 and its mathematical representation, Proc. Natl. Acad. Sci. USA, .(6):275–288, 1920, https:

19. P.A. Stephens, W.J. Sutherland, Consequences of the Allee effect for behaviour, ecology and conservation, Trends Ecol. Evol., 14(10):401–405, 1999, https://doi.org/10.1016/ S0169-5347(99)01684-5.

20. P. Turchin, Complex Population Dynamics a Theoretical/Empirical Synthesis, Princeton Univ. Press, Princeton, NJ, 2003.

21. J. Wei, S. Ruan, Stability and bifurcation in a neurall network model with two delays, Physica D, 130(3–4):255–272, 1999, https://doi.org/10.1016/S0167- 2789(99)00009-3.

22. J. Wei Y. Song, Y. Peng, Bifurcations for a predator–prey system with two delays, J. Math. Anal. Appl., 337(1):466–479, 2008, https://doi.org/10.1016/j.jmaa.2007.04. 001.

23. R. Yang, C. Zhang, Dynamics in a diffusive predator–prey system with a constant prey refuge and delay, Nonlinear Anal., Real World Appl., 31:1–22, 2016, https://doi.org/10. 1016/j.nonrwa.2016.01.005.

24. H. Zhao, X. Zhang, X. Huang, Hopf bifurcation and spatial patterns of a delayed biological economic system with diffusion, Appl. Math. Comput., 266(1):462–480, 2015, https:// doi.org/10.1016/j.amc.2015.05.089.