Hopf bifurcations in a three-species food chain system with multiple delays

1 Introduction

The study on the dynamics of predator-prey system is one of the dominant subjects in ecology and mathematical ecology due to its universal existence and importance. As to our knowledge, delay differential equations exhibit much more complicated dynamics than ordinary differential equations since a time delay could cause a stable equilibrium to become unstable and cause the populations to fluctuate. Thus, time delays of one type or another have been incorporated into mathematical models of population dynamics due to maturation time, capturing time or other reasons. In the last decades, many authors have explored the dynamics of systems with time delay and many interesting results have been obtained [1-13].

However, there may be more species in a habitat and they can construct a food chain. Therefore, it is more realistic to consider a multiple-species predator-prey system. Recently, Baek and Lee [14] proposed the following three-species food chain model with the Lotka-Volterra functional response

where x(t), y(t), z(t) denote the population densities of the lowest-level prey, mid-level predator and top-level predator at time t, respectively. The constant a > 0 is called the intrinsic growth rate of the prey species; b > 0 measures the intraspecific competition of the prey; c₁ > 0 and e₂ > 0 represent the conversion rates of the lowest-level prey to the mid-level predator and the mid-level predator to the top-level predator, respectively; c₁ > 0 and e₁ > 0 measure the hunting of the mid-level predator to the lowest-level prey and the top-level predator to the mid-level predator; d₁ > 0 and d₂ > 0 denote the death rates of the mid-level and top-level predator, respectively.

In system (1), the gestation periods and maturation time of some species are all omitted. However, some species need to take time to have the ability to reproduce and capture food. Based on this fact, by incorporating time delays into the maturation time, Cui and Yan considered the following delayed differential system in [15]

where τ₁≥0 denotes the time from birth to having the ability to predate for top-level predator and τ₂ ≥ 0 represents the maturation time that the mid-level predator can be served as food for the top-level predator, respectively. Delayed models similar to (2) have been investigated widely by many authors and lots of interesting results have been obtained. For more details see [16-18].

In system (2), the author just considered the time for the top-level predator to have the ability to predate and the maturation time for the mid-level predator to be served as food for the top-level predator. But the time for the mid-level predator to have the ability to predate and the maturation time for the lowest-level prey to be served as food for the mid-level predator are omitted. Based on this fact, by incorporating delays into the above terms, we consider a more complicated system with multiple delay

where τ₁≥0 and τ₂≥0 denote the time from birth to having the ability to predate for the mid-level predator and the top-level predator. For a special case, we assume that τ₂+τ₃ and τ₁+τ₃ represent the maturation time for the lowest-level prey to be served as food for the mid-level predator and the maturation time for the mid-level predator to be served as food for the top-level predator, respectively. We address the question how the time delays that we incorporated affect the dynamical properties of the system (3). So the aim of this paper is to study the dynamical behaviors of the system (3), for which we investigate the stability and Hopf bifurcation of a three-species food chain system with multiple delays. We would like to mention that the bifurcation in a predator-prey system with a single or multiple delays had been investigated by many researchers [19-22]. However, to the best of our knowledge, few results for system (3) have been obtained. Therefore, the research of this case is worth considering.

2 Stability of positive equilibrium and existence of local Hopf bifurcations

For convenience, we introduce new variables x₁(t) = x(t), y₁(t) = y(t−τ₁), z₁(t) = z(t−τ₁−τ₂) and assume that τ = τ₁+τ₂+τ₃, so that system (3) can be written as the following system with a single delay:

It is easy to see that the system (3) has a unique positive equilibrium E∗:(x0∗,y0∗,z0∗) provided that the condition (H)ae2c1−d2cc1−d1be2>0

holds, where x0∗=ae2−d2cbe2,y0∗=d2e2,z0∗=ae2c1−d2cc1−d1be2be1e2.

Let u~1(t)=x1(t)−x0∗,u~2(t)=y1(t)−y0∗,u~3(t)=z1(t)−z0∗, and use relations a−bx0∗−cy0∗=0,−d1+c1x0∗−e1z0∗=0,−d2+e2y0∗=0, then system (4) can be rewritten as the following equivalent system

To study the stability of the equilibrium E*, it is sufficient to study the stability of the origin for system (5). The linearized system of system (5) at origin is

The characteristic equation of system (6) is

where A=e1e2y0∗z0∗>0,B=bx0∗>0,C=be1e2x0∗y0∗z0∗>0,D=cc1x0∗y0∗>0.

Next, we will investigate the distribution of roots of Eq.(7). 0bviously, λ = 0 is not a root of Eq.(7). When τ = 0, the characteristic equation becomes

It can be seen that B > 0, B(A+D)−C > 0, and C > 0. Therefore, if follows from the Routh-Hurwitz criteria that all roots of (8) have negative real parts and thus the zero equilibrium of system (5) is asymptotically stable when τ = 0.

Now, we examine when the characteristic equation has pairs of purely imaginary roots. For τ > 0, if iω(ω > 0) is a root of Eq.(7), then ω should satisfy the following equations

Thus

Squaring and adding both the equations of (9), we have

Let z = ω², a1=(AB+BD)2+C2(A+D)2,a2=B2C2−(A+D)4(A+D)2,a3=−2C2,a4=−c4(A+D)2, then Eq. (10) becomes

Thus

where

Let m=8a2−3a1216,n=a13−4a1a2+8a332,D0=n24+m327. Similar to discussion in [24], assume that D₀ > 0, then from the Cardano’s formula for the third-degree algebra equation we know that the equation f(z) = 0 had only one real root z1∗ . If D₀ = 0, then the equation f(z) = 0 has three real roots z₁, z₂ and Z₃ (where z₂ = z₃), and in this case we define z2∗ by max {z₁, z₂}. If D₀ < 0, we know that the equation f(z) = 0 has three different real roots denoted by s1∗,s2∗ and s3∗. In this case,we define z3∗=max{s1∗,s2∗,s3∗}. According to Lemma 2.2 in [24], without loss of generality, we can suppose that Eq. (10) has four positive real roots, denoted by z₁, z₂, z₃, z₄, respectively. Then Eq. (9) should also have four positive real roots ω1=z1,ω2=z2,ω3=z3,ω4=z4. Define

Then (τkj,ωk) are solutions of Eq.(7) and λ = ± i ω_k are a pair of purely imaginary roots of Eq.(7) with τ=τkj. Define τ0=τk00=min1≤k≤4{τk0},ω0=ωk0, where k₀∊ {1, 2, 3, 4}. Then τ⁰ is the first value of τ such that Eq.(7) has purely imaginary roots. In the following discussions, for the sake of convenience, we denoted τkj by τ^j(j = 0,1,2, …) for fixed k ∊ {1,2,3,4}. Let λ(τ) = α(τ)± i ω(τ) be the root of Eq.(7) near τ = τ^j satisfying α(τ^j) = 0, ω(τ^j) = ω₀(j = 0, 1, 2, …).

3 Direction of Hopf bifurcations and stability of bifurcating periodic solutions

In the previous section, we have already obtained that, under certain conditions, the system (4) can undergo Hopf bifurcation at the positive equilibrium E∗:(x0∗,y0∗,z0∗) when τ takes some critical values τ = τ^j, (j = 0, 1, 2, ⋯). In this section, by employing the normal form theory and center manifold theorem introduced by Hassard et al. [23], we shall present the formula determining the direction of the Hopf bifurcation and the stability of bifurcating periodic solutions of (4).

Let u_i(t) = ũ_i(τ t), τ = τ^j + μ, μ ∈ R, then μ = 0 is the Hopf bifurcation value for system (4), dropping the bars for simplification of notations, system (4) becomes the following functional differential equation in C = C([−1,0], R³),

where u(t) = (u₁(t), u₂(t), u₃(t))^T ∈ R³, and L_μ : C → R³, f : R × C → R³ are given by

and

where

By the Riesz representation theorem, there exists a function η(θ, μ) of bounded variation for θ ∈ [− 1, 0], such that

In fact, we can choose

where δ is the Dirac delta function. For ϕ ∈ C([− 1,0], R³), define

and

Then system (12) is equivalent to

where u_t(θ) = u(t + θ) for θ ∈ [− 1, 0]. For ψ ∈ C¹([0,1], (R³)^*), define

and a bilinear inner product

where η(θ) = η(θ, 0). Then A(0) and A^* are adjoint operators. By the discussion in section 2, we know that ±iω₀τ^j are eigenvalues of A(0). Hence, they are also eigenvalues of A^*. We first need to compute the eigenvectors of A(0) and A^* corresponding to iω₀τ_j and −iω₀τ^j, respectively. Suppose q(θ) = (1, q₁, q₂)^Te^iω₀τjθ is the eigenvector of A(0) corresponding to iω₀τ^j, then A(0)q(0) = iω₀τ^jq(0). From the definition of A(0) and (13), (15), (16) we have

For q(− 1) = q(0)e^−iω₀τj, then we obtain

Similarly, we can obtain the eigenvector q∗(s)=D(1,q1∗,q2∗)eiω0τjs of A^* corresponding to −iω₀τ^j, where

In order to assure 〈q^*(s), q(θ)〉 = 1, we need to determine the value of D. By (18), we have

Therefore, we can choose D as

Next we will compute the coordinate to describe the center manifold C₀ at μ = 0. Let u_t be the solution of (17) when μ = 0. Define

On the center manifold C₀, we have

where

z and z¯ are local coordinates for center manifold C₀ in the direction of q^* and q¯∗ . Note that W is real if u_t is real. We only consider real solutions for solutions u_t ∈ C₀ of (17). Since μ = 0, we have

where

Following from (19), (20), (14) and using the algorithms given in [23], we can get the coefficients which will be used to determine the important quantities

where

and E1=(E1(1),E1(2),E1(3))T,E2=(E2(1),E2(2),E2(3))T are both constant vectors.

Similarly to the algorithms given in [25], we can obtain

where

Similarly, we get

where

Thus, we can determine W₂₀(θ) and W₁₁(θ) from (24) and (25). Furthermore, we can compute g₂₁ by (23). Thus we can compute the following values: