Oscillations for a delayed predator–prey model with Hassell–Varley-type functional response

Changjin Xu; Peiluan Li

doi:10.1016/j.crvi.2015.01.002

Population biology/Biologie des populations

Oscillations for a delayed predator–prey model with Hassell–Varley-type functional response

Changjin Xu ¹ ; Peiluan Li ²

¹ Guizhou Key Laboratory of Economics System Simulation, School of Mathematics and Statistics, Guizhou College of Finance and Economics, Guiyang 550004, China
² School of Mathematics and Statistics, Henan University of Science and Technology, Luoyang, Henan 471023, China

Comptes Rendus. Biologies, Volume 338 (2015) no. 4, pp. 227-240.

Abstract

In this paper, a delayed predator–prey model with Hassell–Varley-type functional response is investigated. By choosing the delay as a bifurcation parameter and analyzing the locations on the complex plane of the roots of the associated characteristic equation, the existence of a bifurcation parameter point is determined. It is found that a Hopf bifurcation occurs when the parameter τ passes through a series of critical values. The direction and the stability of Hopf bifurcation periodic solutions are determined by using the normal form theory and the center manifold theorem due to Faria and Maglhalaes (1995). In addition, using a global Hopf bifurcation result of Wu (1998) for functional differential equations, we show the global existence of periodic solutions. Some numerical simulations are presented to substantiate the analytical results. Finally, some biological explanations and the main conclusions are included.

Metadata

Received: 2012-03-01
Accepted: 2015-01-07
Published online: 2015-03-30

PMID

DOI: 10.1016/j.crvi.2015.01.002

Keywords: Predator-prey model, Hassell–Varley functional response, Stability, Hopf bifurcation

Author's affiliations:

Changjin Xu ¹; Peiluan Li ²

¹ Guizhou Key Laboratory of Economics System Simulation, School of Mathematics and Statistics, Guizhou College of Finance and Economics, Guiyang 550004, China
² School of Mathematics and Statistics, Henan University of Science and Technology, Luoyang, Henan 471023, China

@article{CRBIOL_2015__338_4_227_0,
     author = {Changjin Xu and Peiluan Li},
     title = {Oscillations for a delayed predator{\textendash}prey model with {Hassell{\textendash}Varley-type} functional response},
     journal = {Comptes Rendus. Biologies},
     pages = {227--240},
     publisher = {Elsevier},
     volume = {338},
     number = {4},
     year = {2015},
     doi = {10.1016/j.crvi.2015.01.002},
     language = {en},
}

TY  - JOUR
AU  - Changjin Xu
AU  - Peiluan Li
TI  - Oscillations for a delayed predator–prey model with Hassell–Varley-type functional response
JO  - Comptes Rendus. Biologies
PY  - 2015
SP  - 227
EP  - 240
VL  - 338
IS  - 4
PB  - Elsevier
DO  - 10.1016/j.crvi.2015.01.002
LA  - en
ID  - CRBIOL_2015__338_4_227_0
ER  -

%0 Journal Article
%A Changjin Xu
%A Peiluan Li
%T Oscillations for a delayed predator–prey model with Hassell–Varley-type functional response
%J Comptes Rendus. Biologies
%D 2015
%P 227-240
%V 338
%N 4
%I Elsevier
%R 10.1016/j.crvi.2015.01.002
%G en
%F CRBIOL_2015__338_4_227_0

Changjin Xu; Peiluan Li. Oscillations for a delayed predator–prey model with Hassell–Varley-type functional response. Comptes Rendus. Biologies, Volume 338 (2015) no. 4, pp. 227-240. doi : 10.1016/j.crvi.2015.01.002. https://comptes-rendus.academie-sciences.fr/biologies/articles/10.1016/j.crvi.2015.01.002/

Version originale du texte intégral

1 Introduction

Since the pioneering work of Volterra [1] and Lotka [2] in the mid-1920s, there has been increasing interest in investigating the dynamical behaviors of predator–prey models in both ecology and mathematical ecology [3–11]. In particular, one of the important dynamical predator–prey behaviors, such as periodic phenomena and bifurcation has become even more interesting [6–10,12–24]. In 1980, Freeman [25] proposed a most popular predator–prey model with Michaelis–Menten-type functional response:

\{\begin{array}{c} \frac{d x_{1}}{d t} = r x_{1} (1 - \frac{x_{1}}{K}) - \frac{c x_{1} x_{2}}{m + x_{1}}, \\ \frac{d x_{2}}{d t} = x_{2} (- d + \frac{f x_{1}}{m + x_{1}}), \\ x (0) > 0, y (0) > 0, \end{array}

(1)

where x₁, x₂ denote the population of preys and predators at time t, respectively. r, K, c, m, d, and f are positive constants that denote the prey's intrinsic growth rate, carrying capacity, capturing rate, half-saturation constant, predator death rate, maximal predator growth rate, respectively. For more details about the model, the reader is referred to [25].

Considering that in many situations, predators must search and share or compete for food, Arditi and Ginzburg [3] introduced and studied the following ratio-dependent-type functional response model:

\{\begin{array}{c} \frac{d x_{1}}{d t} = r x_{1} (1 - \frac{x_{1}}{K}) - \frac{c x_{1} x_{2}}{m x_{2} + x_{1}}, \\ \frac{d x_{2}}{d t} = x_{2} (- d + \frac{f x_{1}}{m x_{2} + x_{1}}), \\ x (0) > 0, y (0) > 0 \end{array}

(1.2)

Since the functional response depends on the predator density in a different way, Hassel and Varley [26] reconstructed the predator–prey model with Hassell–Varly-type functional response, which takes the following form:

\{\begin{array}{c} \frac{d x_{1}}{d t} = r x_{1} [1 - \frac{x_{1}}{K}] - \frac{c x_{1} x_{2}}{m x_{2}^{γ} + x_{1}}, \\ \frac{d x_{2}}{d t} = x_{2} [- d + \frac{f x_{1}}{m x_{2}^{γ} + x_{1}}], \\ x (0) > 0, y (0) > 0, \end{array}

(1.3)

where γ ∈ (0,1) is called the HV constant. Generally, the consumptions of prey by the predator throughout its past history governs the present birth rate of the predator. Motivated by this point of view, Wang [27] introduced and investigated the periodic solutions to the following delayed predator–prey model:

\{\begin{array}{c} \frac{d x_{1} (t)}{d t} = x_{1} (t) [a (t) - b (t) x_{1} (t - τ (t)) - \frac{c (t) x_{2} (t)}{m x_{2}^{γ} (t) + x_{1} (t)}], \\ \frac{d x_{2} (t)}{d t} = x_{2} (t) [- d (t) + \frac{r (t) x_{1} (t)}{m x_{2}^{γ} + x_{1}}], \\ x (0) 0, y (0) > 0 \end{array}

(1.4)

with the following initial condition:

\{\begin{array}{c} x_{1} (t) = φ (θ), θ \in [- δ, 0], φ (0) = φ_{0} > 0, \\ x_{2} (t) = ψ (θ), θ \in [- δ, 0], ψ (0) = ψ_{0} > 0, \end{array}

(1.5)

where

δ = \sup_{t \in [0, ω]} \{τ (t)\}, φ, ψ \in C ([- δ, 0])

with the norm

∥x∥ = \sup_{t \in [- δ, 0]} |x (t)|

. It is worth pointing out that during the course of the predator–prey interaction when predators do not form groups, one can assume that the HV constant is equal to 1, that is, γ = 1. Moreover, it is more reasonable to incorporate the delay into Hassell–Varly-type functional response. From the point of view of biology, we will consider the following model with delayed Hassell–Varly-type functional response:

\{\begin{array}{c} \frac{d x_{1}}{d t} = x_{1} [a - b x_{1} (t - τ) - \frac{c x_{2} (t - τ)}{m x_{2} (t - τ) + x_{1} (t - τ)}], \\ \frac{d x_{2}}{d t} = x_{2} [- d + \frac{r x_{1} (t - τ)}{m x_{2} (t - τ) + x_{1} (t - τ)}] \end{array}

(1.6)

In this paper, we will devote our attention to investigating the properties of a Hopf bifurcation of system (1.6), that is to say, we shall take the delay τ as the bifurcation parameter and show that when τ passes through a certain critical value, the positive equilibrium loses its stability and a Hopf bifurcation will take place. Furthermore, when the delay τ takes a sequence of critical values containing the above critical value, the positive equilibrium of system (1.6) will undergo a Hopf bifurcation. In particular, by using the normal form theory and the center manifold reduction due to Faria and Maglhalaes [28], the formulae for determining the direction of Hopf bifurcations and the stability of bifurcating periodic solutions are obtained. In addition, the existence of periodic solutions for τ far away from the Hopf bifurcation values is also established by means of the global Hopf bifurcation result of Wu [29].

In order to obtain the main results of our paper, throughout this paper, we assume that the coefficients of system (1.6) satisfy the following condition:

H₁

am² + cd − cr > 0, r > d

This paper is organized as follows. In Section 2, the stability of the positive equilibrium and the existence of a Hopf bifurcation at the positive equilibrium are studied. In Section 3, the direction of Hopf bifurcation and the stability of bifurcating periodic solutions on the center manifold are determined. In Section 4, numerical simulations are carried out to illustrate the validity of the main results. In Section 5, some conditions that guarantee the global existence of the bifurcating periodic solutions to the model are given. Biological explanations and some main conclusions are drawn in Section 6.

2 Stability of the equilibrium and existence of the local Hopf bifurcation

In the section, by analyzing the characteristic equation of the linearized system of system (1.6) at the positive equilibrium, we investigate the stability of the positive equilibrium and the existence of the local Hopf bifurcations occurring at the positive equilibrium.

Considering the biological meaning, we only study the property of a unique positive equilibrium (i.e., coexistence equilibrium). It is easy to see that under the hypothesis (H₁), system (1.6) has a unique positive equilibrium $E_{*} (x_{1}^{*}, x_{2}^{*})$ , where

x_{1}^{*} = \frac{a m^{2} + c d - c r}{a b m}, x_{2}^{*} = \frac{(r - d) (a m^{2} + c d - c r)}{a b d m^{2}}

Let $u_{1} (t) = x_{1} (t) - x_{1}^{*}, u_{2} (t) = x_{2} (t) - x_{2}^{*},$ then, system (1.6) takes the following form:

\{\begin{array}{c} \begin{array}{l} \frac{d u_{1}}{d t} = m_{1} u_{1} (t - τ) + m_{2} u_{2} (t - τ) + m_{3} u_{1}^{2} (t - τ) + m_{4} u_{2}^{2} (t - τ) \\ + m_{5} u_{1} (t - τ) u_{2} (t - τ) + m_{6} u_{1} (t) u_{1} (t - τ) + m_{7} u_{1} (t) u_{2} (t - τ) \\ + m_{8} u_{1} (t) u_{1}^{2} (t - τ) + m_{9} u_{1} (t) u_{2}^{2} (t - τ) + h . o . t ., \end{array} \\ \begin{array}{l} \frac{d u_{2}}{d t} = n_{1} u_{1} (t - τ) + n_{2} u_{2} (t - τ) + n_{3} u_{1}^{2} (t - τ) + n_{4} u_{2}^{2} (t - τ) \\ + n_{5} u_{1} (t - τ) u_{2} (t - τ) + n_{6} u_{1} (t - τ) u_{2} (t) + n_{7} u_{2} (t) u_{2} (t - τ) \\ + n_{8} u_{1} (t - τ) u_{2}^{2} (t - τ) + n_{9} u_{1}^{2} (t - τ) u_{2} (t - τ) + n_{10} u_{2} (t) u_{2}^{2} (t - τ) \\ + n_{11} u_{1} (t - τ) u_{2} (t) u_{2} (t - τ) + n_{12} u_{1}^{2} (t - τ) u_{2} (t) + h . o . t ., \end{array} \end{array}

(2.1)

where

\begin{array}{l} m_{1} = (\frac{c x_{2}^{*}}{m x_{2}^{*} + x_{1}^{*}} - b) x_{1}^{*}, m_{2} = (\frac{m c x_{2}^{*}}{{(m x_{2}^{*} + x_{1}^{*})}^{2}} - \frac{c}{m x_{2}^{*} + x_{1}^{*}}) x_{1}^{*}, \\ m_{3} = - \frac{c x_{1}^{*} x_{2}^{*}}{{(m x_{2}^{*} + x_{1}^{*})}^{3}}, m_{4} = \frac{m c x_{1}^{*}}{{(m x_{2}^{*} + x_{1}^{*})}^{2}} - \frac{m^{2} c x_{2}^{*}}{{(m x_{2}^{*} + x_{1}^{*})}^{3}}, \\ m_{5} = \frac{c}{{(m x_{2}^{*} + x_{1}^{*})}^{2}} - \frac{2 m c x_{2}^{*}}{{(m x_{2}^{*} + x_{1}^{*})}^{2}}, m_{6} = \frac{c x_{2}^{*}}{m x_{2}^{*} + x_{1}^{*}} - b, \\ m_{7} = \frac{m c x_{2}^{*}}{{(m x_{2}^{*} + x_{1}^{*})}^{2}} - \frac{c}{m x_{2}^{*} + x_{1}^{*}}, m_{8} = - \frac{c x_{2}^{*}}{{(m x_{2}^{*} + x_{1}^{*})}^{3}}, \\ m_{9} = \frac{m c}{{(m x_{2}^{*} + x_{1}^{*})}^{2}} - \frac{m^{2} c x_{2}^{*}}{{(m x_{2}^{*} + x_{1}^{*})}^{3}}, n_{1} = \frac{r x_{2}^{*}}{m x_{2}^{*} + x_{1}^{*}} - \frac{r x_{1}^{*} x_{2}^{*}}{{(m x_{2}^{*} + x_{1}^{*})}^{2}}, \\ n_{2} = - \frac{m r x_{1}^{*} x_{2}^{*}}{{(m x_{2}^{*} + x_{1}^{*})}^{2}}, n_{3} = - \frac{r x_{2}^{*}}{{(m x_{2}^{*} + x_{1}^{*})}^{2}}, n_{4} = - \frac{m^{2} r x_{1}^{*} x_{2}^{*}}{{(m x_{2}^{*} + x_{1}^{*})}^{3}}, \\ n_{5} = \frac{2 m r x_{1}^{*} x_{2}^{*}}{{(m x_{2}^{*} + x_{1}^{*})}^{3}} - \frac{m r x_{2}^{*}}{{(m x_{2}^{*} + x_{1}^{*})}^{2}}, n_{6} = \frac{r}{m x_{2}^{*} + x_{1}^{*}} - \frac{r x_{1}^{*}}{{(m x_{2}^{*} + x_{1}^{*})}^{2}}, \\ n_{7} = - \frac{m r x_{1}^{*}}{{(m x_{2}^{*} + x_{1}^{*})}^{2}}, n_{8} = - \frac{m^{2} r x_{2}^{*}}{{(m x_{2}^{*} + x_{1}^{*})}^{3}}, n_{9} = \frac{2 m r x_{2}^{*}}{{(m x_{2}^{*} + x_{1}^{*})}^{3}}, \\ n_{10} = - \frac{m^{2} r x_{1}^{*}}{{(m x_{2}^{*} + x_{1}^{*})}^{3}}, n_{11} = \frac{2 m r x_{1}^{*}}{{(m x_{2}^{*} + x_{1}^{*})}^{3}} - \frac{m r}{{(m x_{2}^{*} + x_{1}^{*})}^{2}}, n_{12} = - \frac{r}{{(m x_{2}^{*} + x_{1}^{*})}^{2}} \end{array}

Then, we obtain the linearized system of (2.1)

\{\begin{array}{c} \frac{d u_{1}}{d t} = m_{1} u_{1} (t - τ) + m_{2} u_{2} (t - τ), \\ \frac{d u_{2}}{d t} = n_{1} u_{1} (t - τ) + n_{2} u_{2} (t - τ) \end{array}

(2.2)

with characteristic equation:

λ^{2} - (m_{1} + n_{2}) λ e^{- λ τ} + (m_{1} n_{2} - m_{2} n_{1}) e^{- 2 λ τ} = 0

(2.3)

In order to investigate the distribution of roots of the transcendental equation (2.3), the following Lemma stated in [30] is helpful.

Lemma 2.1.[30]For the transcendental equation

\begin{array}{l} P (λ, e^{- λ τ_{1}}, \dots, e^{- λ τ_{m}}) = λ^{n} + p_{1}^{(0)} λ^{n - 1} + \dots + p_{n - 1}^{(0)} λ + p_{n}^{(0)} \\ + [p_{1}^{(1)} λ^{n - 1} + \dots + p_{n - 1}^{(1)} λ + p_{n}^{(1)}] e^{- λ τ_{1}} + \dots \\ + [p_{1}^{(m)} λ^{n - 1} + \dots + p_{n - 1}^{(m)} λ + p_{n}^{(m)}] e^{- λ τ_{m}} = 0, \end{array}

as

(τ_{1}, τ_{2}, τ_{3}, \dots, τ_{m})

vary, the sum of orders of the zeros of

P (λ, e^{- λ τ_{1}}, \dots, e^{- λ τ_{m}})

in the open right half plane can change, and only a zero appears on or crosses the imaginary axis.

Regarding τ as the parameter, we can apply Lemma 2.1 to (2.3), which is a special case of

P (λ, e^{- λ τ_{1}}, \dots, e^{- λ τ_{m}})

Obviously, if m₁n₂ ≠ m₂n₁, then, λ = 0 is not a root of (2.3). For τ = 0 the characteristic equation (2.3) becomes:

λ^{2} - (m_{1} + n_{2}) λ + m_{1} n_{2} - m_{2} n_{1} = 0

(2.4)

It is easy to see that Eq. (2.4) have two negative real roots if the following condition:

H₂

m₁ + n₂ < 0, m₁n₂ − m₂n₁ > 0

holds.

Multiplying e^λτ on both sides of (2.3), it is obvious to obtain:

λ^{2} e^{λ τ} - (m_{1} + n_{2}) λ + (m_{1} n_{2} - m_{2} n_{1}) e^{- λ τ} = 0

(2.5)

For ω₀ > 0, iω₀ is a root of (2.5) if and only if

- ω_{0}^{2} e^{i ω_{0} τ} - (m_{1} + n_{2}) i ω + (m_{1} n_{2} - m_{2} n_{1}) e^{- i ω_{0} τ} = 0

(2.6)

Separating the real and imaginary parts of (2.6), we get:

\{\begin{array}{c} (m_{1} n_{2} - m_{2} n_{1} - ω_{0}^{2}) \cos ω_{0} τ = 0, \\ (ω_{0}^{2} + m_{1} n_{2} - m_{2} n_{1}) \sin ω_{0} τ = (m_{1} + n_{2}) ω_{0} \end{array}

(2.7)

If the condition (H₂) holds, we can easily check that cosω₀τ ≠ 0. Then, it follows from (2.7) that:

m_{1} n_{2} - m_{2} n_{1} = ω_{0}^{2}

which leads to:

ω_{0} = \sqrt{m_{1} n_{2} - m_{2} n_{1}}

(2.8)

From the second equation of (2.7), we can easily obtain:

τ_{k} = \frac{1}{ω_{0}} [\arcsin \frac{(m_{1} + n_{2}) \sqrt{m_{1} n_{2} - m_{2} n_{1}}}{2 (m_{1} n_{2} - m_{2} n_{1})} + 2 k π], k = 0,1,2, \dots .

(2.9)

From (2.7), we know that (2.3) has a simple pair of purely imaginary roots ±iω₀ at $τ_{k} (k = 0,1,2, \dots)$

Let $λ_{k} (τ) = α_{k} (τ) + i ω_{k} (τ)$ be the root of Eq. (2.3) satisfying $α_{k} (τ_{k}) = 0, ω_{k} (τ_{k}) = ω_{0}$

Assume that

H₃

$4 {(m_{1} n_{2} - m_{2} n_{1})}^{2} + (m_{1} + n_{2}) \sqrt{4 {(m_{1} n_{2} - m_{2} n_{1})}^{2} - {(m_{1} + n_{2})}^{2}} > 2 {(m_{1} + n_{2})}^{2}$

Then, the following transversality condition holds.

Lemma 2.2. If ( H₁ ), ( H₂ ) and ( H₃ ) are satisfied, then $\frac{d α_{k} (τ)}{d τ} |τ = τ_{k} > 0$

Proof. Differentiating the equation (2.3) with respect to τ leads to:

{[\frac{d λ}{d τ}]}^{- 1} = \frac{2 λ e^{λ τ} + (m_{1} + n_{2}) τ - 2 (m_{1} n_{2} - m_{2} n_{1}) τ e^{- λ τ}}{2 (m_{1} n_{2} - m_{2} n_{1}) τ e^{- λ τ} λ - (m_{1} + n_{2}) λ}

When τ = τ_k, iω₀ is a purely imaginary root of (2.3). We then easily get:

{[\frac{d α_{k} (τ)}{d τ}]}^{- 1} |_{τ = τ_{k}} = \frac{4 ω_{0}^{4} \cos 2 ω_{0} τ_{k} - 2 ω_{0}^{3} (m_{1} + n_{2}) \cos ω_{0} τ_{k}}{{[2 ω_{0}^{3} \sin ω_{0} τ_{k}]}^{2} + {[2 ω_{0}^{3} \cos ω_{0} τ_{k} - (m_{1} + n_{2}) ω_{k}]}^{2}}

(2.10)

From (2.7), we have:

\sin ω_{0} τ_{k} = \frac{m_{1} + n_{2}}{2 ω_{0}}, \cos ω_{0} τ_{k} = \pm \frac{\sqrt{4 ω_{0}^{2} - {(m_{1} + n_{2})}^{2}}}{2 ω_{0}}

Hence,

\begin{array}{l} {[\frac{d α_{k} (τ)}{d τ}]}^{- 1} |_{τ = τ_{k}} = \frac{ω_{0}^{2} {[4 (m_{1} n_{2} - m_{2} n_{1}) - 2 {(m_{1} + n_{2})}^{2}]}^{2}}{{[2 ω_{0}^{3} \sin ω_{0} τ_{k}]}^{2} + {[2 ω_{0}^{3} \cos ω_{0} τ_{k} - (m_{1} + n_{2}) ω_{k}]}^{2}} \\ \pm \frac{ω_{0}^{2} \sqrt{4 (m_{1} n_{2} - m_{2} n_{1}) - {(m_{1} + n_{2})}^{2}}}{{[2 ω_{0}^{3} \sin ω_{0} τ_{k}]}^{2} + {[2 ω_{0}^{3} \cos ω_{0} τ_{k} - (m_{1} + n_{2}) ω_{k}]}^{2}} \end{array}

Under the condition (H₃), we know that

\frac{d α_{k} (τ)}{d τ} |τ = τ_{k} > 0,

completing the proof. Lemma 2.2 implies that the roots λ_k (τ) of characteristic equation (2.3) near τ_k crosses the imaginary axis from the left to the right as τ continuously varies from a number less than τ_k to one greater than τ_k by Rouché’s theorem [31].

Applying Lemma 2.1, we obtain the following results: Lemma 2.3. If (H₁), (H₂) and (H₃) are satisfied, then

(i) if $τ \in [0, τ_{0})$ , all roots Eq. (2.3) have negative real parts;

(ii) if τ = τ₀, all roots of Eq. (2.3) except ±iω₀ have negative real parts;

(iii) if τ ∈ [τ_k, τ_k+1) for k = 0, 1, 2, ..., Eq. (2.3) have 2 (k + 1) roots with positive real parts.

Spectral properties of Eq. (2.3) immediately lead to the properties of the zero solutions to system (2.2), and equivalently, of the positive equilibrium E_∗ for system (1.6).

Theorem 2.4.Suppose that (H₁), (H₂) and (H₃) are satisfied. Then, the positive equilibrium E_∗of system(1.6)is asymptotically stable when $τ \in [0, τ_{0})$ , and unstable when τ > τ₀. Moreover, at τ = τ_k, k = 0, 1, 2, ..., ± iω₀ are a simple pair imaginary roots of (2.3), and (1.6) undergoes a Hopf bifurcation near E_∗.

3 Direction and stability of the Hopf bifurcation

In the previous section, we have obtained conditions for Hopf bifurcations to occur when τ = τ_k. In this section, we will employ the algorithm of Faria and Maglhalaes [28] to compute explicitly the normal forms of system (1.6) on the center manifold. After that, we will investigate the direction of the Hopf bifurcation and stability of the bifurcating periodic orbits from the positive equilibrium E_∗ of system (1.6) at these critical values of τ_k. We know that Eq. (2.3) has a pair of purely imaginary roots ±iω₀ when τ = τ_k and system (1.6) undergoes a Hopf bifurcation from E_∗. Let μ = τ–τ_k, then μ = 0 is the Hopf bifurcation value of system (1.6).

Throughout this section, we refer to Faria and Maglhalaes [28] for the meaning of the notations involved.

Normalizing the delay τ by the time-scaling t → t/τ, then the system (2.1) can be rewritten as a functional differential equation in $C ([- 1,0], R^{2})$ :

\{\begin{array}{c} \begin{array}{l} \frac{d u_{1}}{d t} = τ [m_{1} u_{1} (t - 1) + m_{2} u_{2} (t - 1) + m_{3} u_{1}^{2} (t - 1) + m_{4} u_{2}^{2} (t - 1) \\ + m_{5} u_{1} (t - 1) u_{2} (t - 1) + m_{6} u_{1} (t) u_{1} (t - 1) + m_{7} u_{1} (t) u_{2} (t - 1) \\ + m_{8} u_{1} (t) u_{1}^{2} (t - 1) + m_{9} u_{1} (t) u_{2}^{2} (t - 1) + h . o . t .], \end{array} \\ \begin{array}{l} \frac{d u_{2}}{d t} = τ [n_{1} u_{1} (t - 1) + n_{2} u_{2} (t - 1) + n_{3} u_{1}^{2} (t - 1) + n_{4} u_{2}^{2} (t - 1) \\ + n_{5} u_{1} (t - 1) u_{2} (t - 1) + n_{6} u_{1} (t - 1) u_{2} (t) + n_{7} u_{2} (t) u_{2} (t - 1) \\ + n_{8} u_{1} (t - 1) u_{2}^{2} (t - 1) + n_{9} u_{1}^{2} (t - 1) u_{2} (t - 1) + n_{10} u_{2} (t) u_{2}^{2} (t - 1) \\ + n_{11} u_{1} (t - 1) u_{2} (t) u_{2} (t - 1) + n_{12} u_{1}^{2} (t - 1) u_{2} (t) + h . o . t .] \end{array} \end{array}

(3.1)

Let u = (u₁(t), u₂(t))^T, then, system (3.1) can be rewritten in the following vector form:

\frac{d u}{d t} = L (τ) (u_{t}) + F (u_{t}, τ),

(3.2)

where

L (τ) (ϕ) = τ (\begin{array}{c} m_{1} ϕ_{1} (- 1) + m_{2} ϕ_{2} (- 1) \\ n_{1} ϕ_{1} (- 1) + n_{2} ϕ_{2} (- 1) \end{array}),

F (ϕ, τ) = τ (\begin{array}{c} F_{1} (ϕ, τ) \\ F_{2} (ϕ, τ) \end{array})

where

ϕ = {(ϕ_{1}, ϕ_{2})}^{T} \in C ([- 1,0]), R^{2}),

and

\begin{array}{l} F_{1} (ϕ, τ) = m_{3} ϕ_{1}^{2} (- 1) + m_{4} ϕ_{2}^{2} (- 1) + m_{5} ϕ_{1} (- 1) u_{2} (- 1) + m_{6} ϕ_{1} (0) ϕ_{1} (- 1) \\ + m_{7} ϕ_{1} (0) ϕ_{2} (- 1) + m_{8} ϕ_{1} (0) ϕ_{1}^{2} (- 1) + m_{9} ϕ_{1} (0) ϕ_{2}^{2} (- 1) + h . o . t ., \\ F_{2} (ϕ, τ) = n_{3} ϕ_{1}^{2} (- 1) + n_{4} ϕ_{2}^{2} (- 1) + n_{5} ϕ_{1} (- 1) ϕ_{2} (- 1) + n_{6} ϕ_{1} (- 1) ϕ_{2} (0) \\ + n_{7} ϕ_{2} (0) ϕ_{2} (- 1) + n_{8} ϕ_{1} (- 1) ϕ_{2}^{2} (- 1) + n_{9} ϕ_{1}^{2} (- 1) ϕ_{2} (- 1) \\ + n_{10} ϕ_{2} (0) ϕ_{2}^{2} (- 1) + n_{11} ϕ_{1} (- 1) ϕ_{2} (0) ϕ_{2} (- 1) + n_{12} ϕ_{1}^{2} (- 1) ϕ_{2} (0) + h . o . t .. \end{array}

Obviously, L(τ) is a continuous linear function mapping $C ([- 1,0]), R^{2})$ into R². By the Riesz representation theorem, there exists a 2 × 2 matrix function $η (θ, τ), - 1 \leq θ \leq 0$ , whose elements are of bounded variation such that

ϕ \mapsto Φ (Ψ, ϕ)

(3.3)

In fact, we can choose

η (θ, τ) = - τ (\begin{array}{c} m_{1} & m_{2} \\ n_{1} & n_{2} \end{array}) δ (θ + 1),

(3.4)

where δ denotes Dirac delta function. Then, (3.3) is satisfied. If ϕ is any a given function in

C ([- 1,0]), R^{2})

and u(ϕ) is the unique solution to the linearized equation

\frac{d u}{d t} = L (τ) u_{t}

of Eq. (3.2) with the initial function ϕ at zero, then the solution operator

T (t) : C ([- 1,0], R^{2}) \to C ([- 1,0]), R^{2})

is defined by

T (t) ϕ = u_{t} (ϕ), t \geq 0

(3.5)

From Lemma 7.1.1 in Hale [32], we know that $T (t), t \geq 0$ , is a strongly continuous semigroup of linear transformation on $[0, \infty)$ and the infinitesimal generator A(τ) of T(t), t ≥ 0 is given by

A (τ) ϕ (θ) = \dot{ϕ} (θ) + X_{0} (θ) [L (τ) (ϕ) - \dot{ϕ} (0)] for ϕ \in C ([- 1,0]), R^{2}),

(3.6)

where

X_{0} (θ) = \{\begin{array}{l} 0, & - 1 \leq θ < 0, \\ I, & θ = 0 \end{array}

(3.7)

For $ψ \in C ([- 1,0]), {(R^{2})}^{*})$ , define

A^{*} ψ (s) = - \dot{ψ} (s) + X_{0} (- s) [\int_{- 1}^{0} ψ (- t) d η (t, τ_{k}) + \dot{ψ} (0)]

(3.8)

and a linear inner product

< ψ, ϕ > = \bar{ψ} (0) ϕ (0) - \int_{- 1}^{0} \int_{ξ = 0}^{θ} ψ^{T} (ξ - θ) d η (θ) ϕ (ξ) d ξ,

(3.9)

where

η (θ) = η (θ, 0)

the

A = A (0)

and A^* are adjoint operators. By the discussions in Section 2, we know that ±iω₀τ_k are eigenvalues of A(0), and they are also eigenvalues of A^* corresponding to iω₀τ_k and –iω₀τ_k, respectively. Let

Λ = \{- i ω_{0} τ_{k}, i ω_{0} τ_{k}\}

and denote by P the invariant space of A(τ_k) associated with Λ, where the dimension of P equals to 2. Now, we can decompose the space

C ≔ C ([- 1,0]), R^{2})

as C = P⊕Q by applying the formal adjoint theory for functional differential equations in [32].

Consider the complex coordinates and still denote $C ([- 1,0]), R^{2})$ as C. Suppose that Φ = (Φ₁, Φ₂) is a basis of P and

Φ_{1} (θ) = e^{i ω_{0} τ_{k} θ} {(1, ξ,)}^{T}, Φ_{2} (θ) = \bar{Φ_{1} (θ)}, - 1 \leq θ \leq 0,

where

ξ = \frac{i ω_{0} e^{i ω_{0} τ_{k}} - m_{1}}{m_{2}} .

Also, the two eigenvectors Ψ₁, Ψ₂ of A^* corresponding to the eigenvalues iω₀τ_k, –iω₀τ_k, respectively, construct a basis Ψ = (Ψ₁, Ψ₂)^T of the adjoint space P^* of P and

Ψ_{1} (θ) = E e^{- i ω_{0} τ_{k} s} {(1, η)}^{T}, Ψ_{2} (θ) = \bar{Ψ_{1} (θ)}, 0 \leq s \leq 1,

where

η = - \frac{i ω_{0} e^{i ω_{0} τ_{k}} + m_{1}}{n_{1}},

E = \frac{1}{1 + \bar{ξ} η + m_{1} e^{i ω_{0} τ_{k}} + m_{2} η e^{i ω_{0} τ_{k}} + m_{2} e^{i ω_{0} τ_{k}} + n_{2} \bar{ξ} η e^{i ω_{0} τ_{k}}} .

Thus, $(Ψ, Φ) = ((Ψ_{j}, Φ_{i}), i, j = 1,2) = I_{2},$ where I₂ is a second-order identical matrix. It is known that $\dot{Φ} = Φ B$ , where

B = (\begin{array}{c} i ω_{0} τ_{k} & 0 \\ 0 & - i ω_{0} τ_{k} \end{array}) .

Take the enlarged phase space C by considering the space $B C ≔ \{ϕ : [- 1,0] \to C^{2} | ϕ is continuous on [- 1,0] {and lim}_{θ \to 0^{-}} ϕ (θ) exists\} .$ The projection $ϕ \mapsto Φ (Ψ, ϕ)$ of C upon P, associated with the decomposition C = P⊕Q, is now replaced by π : BC → P such that $π (ϕ + X_{0} α) = Φ [(Ψ, ϕ) + Ψ (0) α]$ .

Thus, we have the decomposition BC = P⊕Kerπ. Using the decomposition $u_{t} = ϕ x (t) + y_{t}, x (t) \in C^{2}, y_{t} \in K e r π \cap C^{1} = Q^{1}$ , and by Theorem 7.6.1 in [32], we can decompose (3.2) as

\{\begin{cases} \frac{d x}{d t} = B x + ψ (0) F_{0} (Φ x + y, μ), \\ \frac{d y}{d t} = A (\tilde{τ}) |Q^{1} y + (I - π) X_{0} F_{0} (Φ x + y, μ), \end{cases}

(3.10)

where

F_{0} (ϕ, μ) = L (μ) (ϕ) + F (ϕ, τ_{k} + μ)

. In view of the Taylor expansion, we denote, respectively,

Ψ (0) F_{0} (Φ x + y, μ)

and

(I - π) X_{0} F_{0} (Φ x + y, μ)

as:

\{\begin{array}{c} Ψ (0) F_{0} (Φ x + y, μ) = \frac{1}{2} f_{2}^{1} (x, y, μ) + \frac{1}{3} f_{3}^{1} (x, y, μ) + \dots, \\ (I - π) X_{0} F_{0} (Φ x + y, μ) = \frac{1}{2} f_{2}^{2} (x, y, μ) + \frac{1}{3} f_{3}^{2} (x, y, μ) + \dots, \end{array}

(3.11)

where

f_{j}^{1} (x, y, μ)

and

f_{j}^{2} (x, y, μ)

are homogeneous polynomials in (x, y, μ) of degree j, j = 2, 3,…, with coefficients in C² and Kerπ, respectively. The normal form method gives a normal form on the center manifold of the origin for (3.10) as

\dot{x} = B x + \frac{1}{2} g_{2}^{1} (x, 0, μ) + \frac{1}{3} g_{3}^{1} (x, 0, μ) + \dots,

(3.12)

where

g_{j}^{1} (x, 0, μ)

are homogeneous polynomials in (x, μ) of degree j, j = 2, 3,…

In what follows, we first define the operators $M_{j}^{1}$ as

M_{j}^{1} (p) (x, μ) = D_{x} p (x, μ) B x - B p (x, μ), j \geq 2

(3.13)

In particular, $M_{j}^{1} (μ^{l} x^{q} e_{k}) = i ω_{0} τ_{k} (q_{1} - q_{2} + {(- 1)}^{k}) μ^{l} x^{q} e_{k}, l + q_{1} + q_{2} = j, k = 1,2 for j = 2,3, q = (q_{1}, q_{2}) \in N_{0}^{2}, l \in N_{0} and \{e_{1}, e_{2}\}$ is the canonical basis for C². Therefore, we have

\begin{array}{l} Ker (M_{2}^{1}) = span \{(\begin{array}{c} x_{1} μ \\ 0 \end{array}), (\begin{array}{c} 0 \\ x_{2} μ \end{array})\}, \\ Ker (M_{3}^{1}) = span \{(\begin{array}{c} x_{1}^{2} x_{2} \\ 0 \end{array}), (\begin{array}{c} x_{1} μ^{2} \\ 0 \end{array}), (\begin{array}{c} 0 \\ x_{1} x_{2}^{2} \end{array}), (\begin{array}{c} 0 \\ x_{2} μ^{2} \end{array})\} \end{array}

By (3.10), we have

f_{2}^{1} (x, y, μ) = 2 Ψ (0) [L (μ) (Φ x + y) + F_{2} (Φ x + y, τ_{k})]

(3.14)

Noting that L(μ) = μ / τ_kL(τ_k), we have

f_{2}^{1} (x, 0, u) = (\begin{array}{c} 2 A_{1} x_{1} μ + 2 A_{1} μ x_{2} + a_{20} x_{1}^{2} + 2 a_{11} x_{1} x_{2} + a_{02} x_{2}^{2} \\ 2 {\bar{A}}_{1} x_{1} μ + 2 {\bar{A}}_{2} x_{2} μ x_{2} + {\bar{a}}_{20} x_{1}^{2} + 2 {\bar{a}}_{11} x_{1} x_{2} + {\bar{a}}_{02} x_{2}^{2} \end{array})

(3.15)

where

\begin{array}{l} A_{1} = E [m_{1} e^{- i ω_{0} τ_{k}} + m_{2} e^{- i ω_{0} τ_{k}} ξ + η (n_{1} e^{- i ω_{0} τ_{k}} + n_{2} e^{- i ω_{0} τ_{k}} ξ)], \\ A_{2} = E [m_{1} e^{i ω_{0} τ_{k}} + m_{2} e^{i ω_{0} τ_{k}} \bar{ξ} + η (n_{1} e^{i ω_{0} τ_{k}} + n_{2} e^{i ω_{0} τ_{k}} \bar{ξ})], \\ a_{20} = 2 E τ_{k} [m_{3} e^{- 2 i ω_{0} τ_{k}} + m_{4} e^{- 2 i ω_{0} τ_{k}} ξ^{2} + m_{5} e^{- 2 i ω_{0} τ_{k}} ξ + m_{6} e^{- i ω_{0} τ_{k}} + m_{7} e^{- i ω_{0} τ_{k}} ξ], \\ a_{11} = 2 E τ_{k} [m_{3} + m_{4} | ξ |^{2} + m_{5} Re \{ξ\} + m_{6} Re \{e^{i ω_{0} τ_{k}}\} + m_{7} e^{i ω_{0} τ_{k}} \bar{ξ}], \\ a_{02} = 2 E τ_{k} [m_{3} e^{2 i ω_{0} τ_{k}} + m_{4} e^{2 i ω_{0} τ_{k}} + m_{5} e^{2 i ω_{0} τ_{k}} \bar{ξ} + m_{6} e^{i ω_{0} τ_{k}} + m_{7} e^{i ω_{0} τ_{k}} \bar{ξ}] \end{array}

(3.16)

Since the second-order terms in (μ, x) on the center manifold are given by

\frac{1}{2} g_{2}^{1} (x, 0, μ) = \frac{1}{2} {Proj}_{Ker (M_{2}^{1})} f_{2}^{1} (x, 0, μ),

(3.17)

it follows that

\frac{1}{2} g_{2}^{1} (x, 0, μ) = (\begin{array}{c} A_{1} x_{1} μ \\ {\bar{A}}_{1} x_{2} μ \end{array}),

(3.18)

where A₁ defined by (3.16).

In the following, we shall compute the cubic term $g_{3}^{1} (x, 0, μ)$ . First, we note that

g_{3}^{1} (x, 0, μ) \in K e r (M_{3}^{1}) = Span \{(\begin{array}{c} x_{1}^{2} x_{2} \\ 0 \end{array}), (\begin{array}{c} x_{1} μ^{2} \\ 0 \end{array}), (\begin{array}{c} 0 \\ x_{1} x_{2}^{2} \end{array}), (\begin{array}{c} 0 \\ x_{2} μ^{2} \end{array})\}

However, the terms $O (0 |x| μ^{2})$ are irrelevant to determine the generic Hopf bifurcation. Hence, it is only needed to compute the coefficients of $(\begin{array}{c} x_{1}^{2} x_{2} \\ 0 \end{array})$ and $(\begin{array}{c} 0 \\ x_{1} x_{2}^{2} \end{array})$ . Let

s ≔ span \{(\begin{array}{c} x_{1} μ \\ 0 \end{array}), (\begin{array}{c} 0 \\ x_{2} μ \end{array})\},

then we have

\frac{1}{3} g_{3}^{1} (x, 0, μ) = \frac{1}{3!} {Proj}_{K e r (M_{3}^{1})} {\bar{f}}_{3}^{1} (x, 0, μ) = \frac{1}{3!} {Proj}_{s} {\bar{f}}_{3}^{1} (x, 0,0) + O (0 |x| μ^{2}),

where

{\bar{f}}_{3}^{1} (x, 0,0) = f_{3}^{1} (x, 0,0) + \frac{3}{2} {[(D_{x} f_{2}^{1}) U_{2}^{1} - (D_{x} U_{2}^{1} g_{2}^{1})]}_{(x, 0,0)} + \frac{3}{2} {[(D_{y} f_{2}^{1}) h]}_{(x, 0,0)}

is the third-order term of the equation, which is obtained after computing the second-order terms of the normal form,

U_{2}^{1} (x, 0)

is the solution to equation

M_{2}^{1} U_{2}^{1} (x, 0) = f_{2}^{1} (x, 0,0)

, and h = (h₁, h₂)^T is a second homogeneous polynomial in (x₁, x₁, μ) with coefficients in Q¹.

From (3.1) and (3.2) and from (3.18) we can easily see that $g_{2}^{1} (x, 0, 0) = 0$ and

{Proj}_{s} f_{3}^{1} (x, 0,0) = (\begin{array}{c} 3 a_{21} x_{1}^{2} x_{2} \\ 3 {\bar{a}}_{21} x_{1} x_{2}^{2} \end{array}),

where

\begin{array}{l} a_{21} = E (2 m_{8} + 2 m_{9} | ξ |^{2}) + E η (2 n_{8} | ξ |^{2} e^{- i ω_{0} τ_{k}} + 2 n_{9} ξ e^{- i ω_{0} τ_{k}} \\ + 2 n_{10} | ξ |^{2} ξ + n_{11} (| ξ |^{2} e^{- i ω_{0} τ_{k}} + ξ^{2} + \bar{ξ}) + 2 n_{12} ξ \end{array}

Therefore, we have ${[(D_{x} U_{2}^{1}) g_{2}^{1}]}_{(x, 0,0)} = 0$ . In the sequel, we only need to compute $U_{2}^{1} (x, 0)$ and h(x)(θ).

From (3.15), we have:

f_{2}^{1} (x, 0,0) = (\begin{array}{c} a_{20} x_{1}^{2} + 2 a_{11} x_{1} x_{2} + a_{02} x_{2}^{2} \\ {\bar{a}}_{02} x_{1}^{2} + 2 {\bar{a}}_{11} x_{1} x_{2} + {\bar{a}}_{20} x_{2}^{2} \end{array})

In view of the definition of $M_{2}^{1}$ , the equation $M_{2}^{1} U_{2}^{1} (x, 0) = f_{2}^{1} (x, 0,0)$ can be written as the following partial differential equations:

\{\begin{array}{c} x_{1} \frac{\partial u_{1}}{\partial x_{1}} - x_{2} \frac{\partial u_{1}}{\partial x_{2}} - u_{1} = \frac{1}{i ω_{0} τ_{k}} (a_{20} x_{1}^{2} + 2 a_{11} x_{1} x_{2} + a_{02} x_{2}^{2}), \\ x_{1} \frac{\partial u_{2}}{\partial x_{1}} - x_{2} \frac{\partial u_{2}}{\partial x_{2}} + u_{2} = \frac{1}{i ω_{0} τ_{k}} ({\bar{a}}_{02} x_{1}^{2} + 2 {\bar{a}}_{11} x_{1} x_{2} + {\bar{a}}_{20} x_{2}^{2}) \end{array}

(3.19)

From (3.19), we can easily obtain:

U_{2}^{1} (x, 0) = (\begin{array}{c} \frac{1}{i ω_{0} τ_{k}} (a_{20} x_{1}^{2} - 2 a_{11} x_{1} x_{2} - \frac{1}{3} a_{02} x_{2}^{2}) \\ \frac{1}{i ω_{0} τ_{k}} (\frac{1}{3} {\bar{a}}_{02} x_{1}^{2} + 2 {\bar{a}}_{11} x_{1} x_{2} - {\bar{a}}_{20} x_{2}^{2}) \end{array})

Thus, we obtain:

{Proj}_{s} {[[D_{x} f_{2}^{1}) U_{2}^{1}]}_{(x, 0,0)} = (\begin{array}{c} \frac{2 i}{ω_{0} τ_{k}} (a_{20} a_{11} - 2 | a_{11} |^{2} - \frac{1}{3} | a_{02} |^{2}) x_{1}^{2} x_{2}) \\ \frac{- 2 i}{ω_{0} τ_{k}} ({\bar{a}}_{20} {\bar{a}}_{11} - 2 | a_{11} |^{2} - \frac{1}{3} | a_{02} |^{2}) x_{1}^{2} x_{2}) \end{array})

In what follows, we shall compute ${Proj}_{s} {[D_{x} f_{2}^{1}) h]}_{(x, 0,0)}$ . From (3.14), we know

\begin{array}{l} f_{2}^{1} (x, y, 0) = 2 Ψ (0) F (Φ x + y, τ_{k}) \\ = 2 τ_{k} (\begin{array}{c} E (m_{3} e_{1}^{2} + m_{4} e_{2}^{2} + m_{5} e_{1} e_{2} + m_{6} e_{1} e_{3} + m_{7} e_{2} e_{3}) \\ + E η (n_{3} e_{1}^{2} + n_{5} e_{2}^{2} + n_{6} e_{1} e_{4} + n_{7} e_{2} e_{4}) \\ \bar{E} (m_{3} e_{1}^{2} + m_{4} e_{2}^{2} + m_{5} e_{1} e_{2} + m_{6} e_{1} e_{3} + m_{7} e_{2} e_{3}) \\ + \bar{E} \bar{η} (n_{3} e_{1}^{2} + n_{5} e_{2}^{2} + n_{6} e_{1} e_{4} + n_{7} e_{2} e_{4}) \end{array}), \end{array}

(3.20)

where

\begin{array}{l} e_{1} = e^{- i ω_{0} τ_{k}} x_{1} + e^{i ω_{0} τ_{k}} x_{2} + y_{1} (- 1), \\ e_{2} = e^{- i ω_{0} τ_{k}} ξ x_{1} + e^{i ω_{0} τ_{k}} \bar{ξ} x_{2} + y_{2} (- 1), \\ e_{3} = x_{1} + x_{2} + y_{1} (0), \\ e_{4} = x_{1} + x_{2} + y_{2} (0) \end{array}

since h = (h⁽¹⁾, h⁽²⁾)^T is a second-order homogeneous polynomial in (x₁, x₂, μ) with coefficients in Q¹. Hence, we can let:

h = h_{110} x_{1} x_{2} + h_{101} x_{1} μ + h_{011} x_{2} μ + h_{200} x_{1}^{2} + h_{020} x_{2}^{2} + h_{002} μ^{2}

(3.21)

Thus, from (3.20), we obtain

{[D_{y} f_{2}^{1}) h]}_{(x, 0,0,)} = 2 τ_{k} (\begin{array}{c} K_{1} \\ K_{2} \end{array}),

where

\begin{array}{l} K_{1} = E [2 m_{3} e_{1}^{0} h^{(1)} (- 1) + m_{5} e_{2}^{0} h^{(1)} (- 1) + e_{1}^{0} h^{(1)} (0)] \\ + m_{7} e_{2}^{0} h^{(1)} (0)] + E η [2 n_{3} e_{1}^{0} h^{(1)} (- 1) + n_{6} e_{4}^{0} h^{(1)} (- 1)] + E [2 m_{4} e_{2}^{0} h^{(2)} (- 1), \\ + m_{5} e_{1}^{0} h^{(2)} (- 1) + m_{7} e_{3}^{0} h^{(2)} (- 1)] + E η [2 n_{5} e_{2}^{0} h^{(2)} (- 1) + n_{6} e_{1}^{0} h^{(2)} (0) \\ + n_{7} (e_{4}^{0} h^{(2)} (- 1) + e_{2}^{0} h^{(2)} (0))], \\ K_{2} = \bar{E} [2 m_{3} e_{1}^{0} h^{(1)} (- 1) + m_{5} (e_{2}^{0} h^{(1)} (- 1) + e_{1}^{0} h^{(1)} (- 1)) + m_{6} e_{3}^{0} h^{(1)} (- 1) \\ + m_{7} e_{2}^{0} h^{(1)} (0)] + \bar{E} \bar{η} [2 n_{3} e_{1}^{0} h^{(1)} (- 1) + n_{6} e_{4}^{0} h^{(1)} (- 1)] + \bar{E} [2 m_{4} e_{2}^{0} h^{(2)} (- 1) \\ + m_{5} e_{1}^{0} h^{(2)} (- 1) + m_{7} e_{3}^{0} h^{(2)} (- 1)] + \bar{E} \bar{η} [2 n_{5} e_{2}^{0} h^{(2)} (- 1) + n_{6} e_{1}^{0} h (2) (0) \\ + n_{7} (e_{4}^{0} h^{(2)} (- 1) + e_{2}^{0} h^{(2)} (0)], \end{array}

where

\begin{array}{l} e_{1}^{0} = e^{- i ω_{0} τ_{k}} x_{1} + e^{i ω_{0} τ_{k}} x_{2}, e_{2}^{0} = e^{- i ω_{0} τ_{k}} ξ x_{1} + e^{i ω_{0} τ_{k}} \bar{ξ} x_{2}, \\ e_{3}^{0} = x_{1} + x_{2}, e_{4}^{0} = x_{1} + x_{2} \end{array}

Therefore,

{[D_{y} f_{2}^{1}) h]}_{(x, 0,0,)} = (\begin{array}{c} 2 B_{3} x_{1}^{2} x_{2} \\ 2 {\bar{B}}_{3} x_{1} x_{2}^{2} \end{array}),

where

B_{3} = - m E τ_{k} (ξ_{3} h_{110}^{(1)} (0) + {\bar{ξ}}_{3} h_{200}^{(1)} (0) + h_{110}^{(3)} (0) + h_{200}^{(3)} (0))

Since h₁₁₀(θ) and h₂₀₀(θ) for $θ \in [- 1,0]$ appear in B₃, we still need to compute them.

Following [28], we know that h = h (θ; x₁, x₂, μ) is the unique solution in the linear space of homogeneous polynomials of degree 2 in 3 real variable (x, μ) = (x₁, x₂, μ) of the equation

(M_{2}^{2} (h) (x, μ) = 2 (I - π) X_{0} [L (μ) (Φ x) + F (Φ x, τ_{k})],

(3.22)

since

(M_{2}^{2} (h) (x, μ) = D_{x} h (x, μ) B x - A (τ_{k}) |Q^{1} (h (x, μ))

(3.23)

Combining the definition (3.6) of the operator A(τ), we can obtain:

\begin{array}{l} D_{x} h (x, μ) B x - \dot{h} (x, μ) - X_{0} (θ) [L (τ) (h (x, μ)) - \dot{h} (x, μ) (0)] \\ = 2 (I - π) X_{0} [L (μ) (Φ x) + F (Φ x, τ_{k})] \end{array}

For the sake of simplicity, let:

h = h_{110} (θ) x_{1} x_{2} + h_{200} (θ) x_{1}^{2} + h_{020} (θ) x_{2}^{2}

(3.24)

Then, h₀ can be evaluated by the system

{\dot{h}}_{0} (x) - D_{x} h_{0} (x, μ) B x = 2 Φ Ψ (0) [L (0) (Φ x) + F (Φ x, τ_{k})],

(3.25)

{\dot{h}}_{0} (x) - L (τ_{k}) (h_{0} (x)) = 2 [L (0) (Φ x) + F (Φ x, τ_{k})],

(3.26)

where

{\dot{h}}_{0}

denote the derivative of h₀ with respect to θ.

In view of (3.14), (3.15), (3.25) and (3.26), we know that $h_{110} = {(h_{110}^{1}, h_{110}^{2})}^{T}$ and $h_{200} = {(h_{200}^{1}, h_{200}^{2})}^{T}$ are the solution to the following two equations:

\{\begin{array}{c} {\dot{h}}_{110} = 2 (a_{11} Φ_{1} + {\bar{a}}_{11} Φ_{2}), \\ {\dot{h}}_{110} (0) - L (τ_{k}) (h_{0} (x)) = 4 τ_{k} (\begin{array}{c} p_{1} \\ p_{2} \end{array}), \end{array}

(3.27)

\{\begin{array}{c} {\dot{h}}_{200} - 2 i ω_{0} τ_{k} h_{200} = a_{20} Φ_{1} + {\bar{a}}_{02} Φ_{2}), \\ {\dot{h}}_{200} (0) - L (τ_{k}) (h_{200}) = τ_{k} (\begin{array}{c} q_{1} \\ q_{2} \end{array}), \end{array}

(3.28)

respectively, where

\begin{array}{l} p_{1} = m_{3} + m_{4} {|ξ|}^{2} + m_{5} Re \{ξ\} + m_{6} Re \{e^{i ω_{0} τ_{k}}\} + m_{7} Re \{ξ e^{- i ω_{0} τ_{k}}\} \\ p_{2} = n_{3} + n_{4} {|ξ|}^{2} + n_{5} Re \{ξ\} + n_{6} Re \{ξ e^{- i ω_{0} τ_{k}}\} + n_{7} Re \{ξ e^{- i ω_{0} τ_{k}}\}, \\ q_{1} = m_{3} e^{- 2 i ω_{0} τ_{k}} + m_{4} e^{- 2 i ω_{0} τ_{k}} ξ^{2} + m_{5} ξ + m_{6} e^{- i ω_{0} τ_{k}} + m_{7} ξ e^{- i ω_{0} τ_{k}}, \\ q_{2} = n_{3} e^{- 2 i ω_{0} τ_{k}} + n_{4} e^{- 2 i ω_{0} τ_{k}} ξ^{2} + n_{5} ξ + n_{6} ξ e^{- i ω_{0} τ_{k}} + n_{7} ξ e^{- i ω_{0} τ_{k}} \end{array}

Solving the Eq. (3.27) and (3.28), we have

\{\begin{array}{c} h_{110} (θ) = \frac{2}{i ω_{0} τ_{k}} (a_{11} Φ_{1} - {\bar{a}}_{11} Φ_{2}) + C, \\ h_{200} (θ) = - \frac{a_{20}}{i ω_{0} τ_{k}} Φ_{1} - \frac{{\bar{a}}_{20}}{3 i ω_{0} τ_{k}} Φ_{2} + D e^{2 i ω_{0} τ_{k}}, \end{array}

(3.29)

where C = (c₁, c₂)^T, D = (d₁, d₂)^T are the solution to the following linear algebra equations:

(\begin{array}{c} m_{1} e^{- i ω_{0} τ_{k}} & m_{2} e^{- i ω_{0} τ_{k}} \\ n_{1} e^{- i ω_{0} τ_{k}} & n_{2} e^{- i ω_{0} τ_{k}} \end{array}) (\begin{array}{c} c_{1} \\ c_{2} \end{array}) = - 4 (\begin{array}{c} p_{1} \\ p_{2} \end{array}),

(3.30)

(\begin{array}{c} 2 i ω_{0} - m_{1} e^{2 i ω_{0} τ_{k}} & - m_{2} e^{2 i ω_{0} τ_{k}} \\ - n_{1} e^{2 i ω_{0} τ_{k}} & 2 i ω_{0} - m_{2} e^{2 i ω_{0} τ_{k}} \end{array}) (\begin{array}{c} d_{1} \\ d_{2} \end{array}) = (\begin{array}{c} q_{1} \\ q_{2} \end{array}),

(3.31)

respectively. We know that

\frac{1}{3} g_{3}^{1} (x, 0,0) = (\begin{array}{c} A_{3} x_{1}^{2} x_{2} \\ {\bar{A}}_{3} x_{1} x_{2}^{2} \end{array}),

(3.32)

where

A_{3} = \frac{i}{2 ω_{0} τ_{k}} (a_{20} a_{11} - 2 {|a_{11}|}^{2} - \frac{1}{3} {|a_{02}|}^{2}) - \frac{1}{2} B_{3}

(3.33)

Accordingly, the normal form (3.12) of (3.10) has the following form

= B x + (\begin{array}{c} A_{1} x_{1} μ \\ {\bar{A}}_{1} x_{2} μ \end{array}) + (\begin{array}{c} A_{3} x_{1}^{2} x_{2} \\ {\bar{A}}_{3} x_{1} x_{2}^{2} \end{array}) + O (|x| μ^{2} + {|x|}^{4})

The normal form (3.12) relative to P can be written in real coordinates (ω₁, ω₂) through the change of variables x₁ = ω₁ − iω₂, x₂ = ω₁ + iω₂. Setting ω₁ = ρ cosν, ω₂ = ρ sinν, then this form becomes

\{\begin{array}{c} \dot{ρ} = k_{1} μ ρ + k_{2} ρ^{3} + O (μ^{2} ρ + {|(ρ, μ)|}^{4}), \\ \dot{ν} = - σ_{k} + O (|(ρ, μ)|), \end{array}

(3.34)

where

k_{1} = Re \{A_{1}\}, k_{2} = Re \{A_{3}\}

. Following [33], we know that the sign of k₁k₂ determines the direction of the Hopf bifurcation and the sign of k₂ determines the stability of the nontrivial periodic solution bifurcating from the Hopf bifurcation. The Hopf bifurcation is supercritical when k₁k₂ < 0 and subcritical if k₁k₂ > 0. In addition, the bifurcating periodic solution on the center manifold is stable provided that k₂ < 0 and unstable if k₂ > 0.

Summarizing the above analysis, we have the following result.

Theorem 3.1.The flow of Eq.(3.2)with μ = 0 on the center manifold of the origin is given by (3.34). Hopf bifurcation is supercritical if k₁k₂ < 0 and subcritical if k₁k₂ > 0. Moreover, the nontrivial periodic solution is stable if k₂ < 0 and unstable if k₂ > 0.

4 Numerical examples

In this section, we present some numerical results for some particular values of the parameters associated with the model system (1.6). We consider the system (1.6) with a = 1, b = 2, c = 0.3, m = 0.5, d = 0.8, r = 2. That is,

\{\begin{array}{c} \frac{d x_{1}}{d t} = x_{1} [1 - 2 x_{1} (t - τ) - \frac{0.3 x_{2} (t - τ)}{0.5 x_{2} (t - τ) + x_{1} (t - τ)}], \\ \frac{d x_{2}}{d t} = x_{2} [- 0.8 + \frac{2 x_{1} (t - τ)}{0.5 x_{2} (t - τ) + x_{1} (t - τ)}], \end{array}

(4.35)

which has a positive equilibrium

E_{*} = (0.32,0.96)

. By means of software Matlab, we obtain:

\begin{array}{l} ω_{0} = 0.9254, τ_{0} = 1.940, ξ = 10.2032 - 5.3126 i, η = 0.1073 + 0.2435 i, \\ E = 0.0107 - 0.04131 i, A_{1} = 0.1912 + 0.3103 i, A_{2} = - 0.2956 - 0.2391 i, \\ a_{20} = 0.7443 - 0.3476 i, a_{11} = 0.0252 - 0.0395 i, a_{02} = - 0.7005 - 0.4669 i, \\ c_{1} = 1.0411 + 0.0012 i, c_{2} = - 0.3589 - 0.8993 i, d_{1} = 2.1251 + 0.9696 i, \\ d_{2} = 0.4998 + 0.5933 i, h_{110}^{(1)} (0) = 0.2013 - 10.0882 i, h_{200}^{(1)} (0) = 2.3266 + 7.0852 i, \\ B_{3} = 0.3092 - 0.5112 i, A_{3} = - 0.1377 - 1.3126 i \end{array}

Then, the stability determining quantities for Hopf bifurcating periodic solutions are given by k₁ = 0.1912, k₂ = −0.1377. From Theorem 2.4, we know that the positive equilibrium $E_{*} (0.32,0.96)$ is asymptotically stable when $τ \in [0, 1.94)$ and is unstable when τ > 1.94. The numerical simulations for τ = 1.9 and τ = 2 are shown on Figs. 1–8, respectively. From Theorem 3.1, we know that system (4.35) with τ = 1.94 has a supercritical Hopf bifurcation and the nontrivial periodic solution bifurcating from the Hopf bifurcation of (4.35) with τ = 1.94 is orbitally asymptotically stable on the center. Moreover, all the roots of Eq. (2.3) with τ = 1.94, except ±0.9254i, have negative real parts. Thus, the center manifold theory implies that the stability of the periodic solutions projected on the center manifold coincides with the stability of the periodic solutions in the whole phase space.

Fig. 1
Trajectory portrait and phase portrait of system (4.35) with τ = 1.9 < τ₀ ≈ 1.94. The positive equilibrium $E_{*} (0.32,0.96)$ is asymptotically stable. The initial value is (0.2, 1.5).

Fig. 2
Trajectory portrait and phase portrait of system (4.35) with τ = 1.9 < τ₀ ≈ 1.94. The positive equilibrium $E_{*} (0.32,0.96)$ is asymptotically stable. The initial value is (0.2, 1.5).

Fig. 3
Trajectory portrait and phase portrait of system (4.35) with τ = 1.9 < τ₀ ≈ 1.94. The positive equilibrium $E_{*} (0.32,0.96)$ is asymptotically stable. The initial value is (0.2, 1.5).

Fig. 4
Trajectory portrait and phase portrait of system (4.35) with τ = 1.9 < τ₀ ≈ 1.94. The positive equilibrium $E_{*} (0.32,0.96)$ is asymptotically stable. The initial value is (0.2, 1.5).

Fig. 5
Trajectory portrait and phase portrait of system (4.35) with τ = 2 < τ₀ ≈ 1.94. Hopf bifurcation occurs from the positive equilibrium $E_{*} (0.32,0.96)$ . The initial value is (0.2, 1.5).

Fig. 6
Trajectory portrait and phase portrait of system (4.35) with τ = 2 < τ₀ ≈ 1.94. Hopf bifurcation occurs from the positive equilibrium $E_{*} (0.32,0.96)$ . The initial value is (0.2, 1.5).

Fig. 7
Trajectory portrait and phase portrait of system (4.35) with τ = 2 < τ₀ ≈ 1.94. Hopf bifurcation occurs from the positive equilibrium $E_{*} (0.32,0.96)$ . The initial value is (0.2, 1.5).

Fig. 8
Trajectory portrait and phase portrait of system (4.35) with τ = 2 < τ₀ ≈ 1.94. Hopf bifurcation occurs from the positive equilibrium $E_{*} (0.32,0.96)$ . The initial value is (0.2, 1.5).

5 Global continuation of local Hopf bifurcations

In the earlier sections we have established that the model system (1.6) undergoes a Hopf bifurcation at $E_{*} (x_{1}^{*}, x_{2}^{*})$ for τ = τ_k and also investigated the stability of bifurcating periodic solutions. We all know that periodic solutions can arise through the Hopf bifurcation in delay differential equations. However, these bifurcating periodic solutions are generally local, i.e., they only exist in a small neighborhood of the critical value. Therefore, we want to know that whether these nonconstant periodic solutions obtained through local Hopf bifurcations exist globally or not. In this section, we will consider the global continuation of periodic solutions bifurcating from the positive equilibrium $E_{*} (x_{1}^{*}, x_{2}^{*})$ of system (1.6).

Throughout this section, we closely follow the notations in Wu [29]. For the simplification of notations, setting $z (t) = {(z_{1} (t), z_{2} (t))}^{T} = {(x_{1} (t), x_{2} (t))}^{T}$ , we can rewrite system (1.6) as the following form

\dot{z} (t) = F (z_{t}, τ, p),

(5.1)

where

z_{t} (θ) = {(z_{1 t} (θ), z_{2 t} (θ))}^{T} = (z_{1} (t + θ), z_{2} (t + θ)) \in ([- τ, 0], R^{2})

. It is easy to see that if the condition (H₁) holds, system (5.1) has a positive equilibrium

E_{*} (x_{1}^{*}, x_{2}^{*})

.

Following the work of Wu [29], we need to define:

\begin{array}{l} X = C ([- τ, 0], R^{2}), \\ Γ = C l \{(z, τ, p) \in X \times R \times R^{+}; z is a nonconstant periodic solution to (5 .1)\}, \\ N = \{(\bar{z}, τ, p); F (\bar{z}, τ, p) = 0\} . \end{array}

Let $\bar{z}$ be the equilibrium of system (5.1). Then, the characteristic matrix of (5.1) at the equilibrium $\bar{z}$ takes the following form

Δ (\bar{z}, τ, p) (λ) = λ I d - D F (\bar{z}, τ, p) (e^{λ} \cdot I d),

where Id is the identity matrix and

D F (\bar{z}, τ, p)

is the Fréchet derivative of F with respect to z_t evaluated at

(\bar{z}, τ, p)

. From Wu [29], we know that

(\bar{z}, τ, p)

is called a center if

(\bar{z}, τ, p) \in N

and

Δ (\bar{z}, τ, p) (λ) = 0

. A center

(\bar{z}, τ, p)

is said to be isolated if it is the only center in some neighborhood of

(\bar{z}, τ, p)

.

For the benefit of the readers, we first state the global Hopf bifurcation theory due to Wu [29] for functional differential equations.

Lemma 5.1. Assume that $(\bar{z}, τ, p)$ is an isolated center satisfying the hypotheses (A ₁ ) – (A ₄ ) in Wu [7] . Denote by

$C (\bar{z}, τ, p)$ the connected component of $(\bar{z}, τ, p)$ in Γ. Then, either

(i) $C (\bar{z}, τ, p)$ is unbounded, or

(ii) $C (\bar{z}, τ, p)$ is bounded, $C (\bar{z}, τ, p) \cap Γ$ is finite and

\sum_{(\bar{z}, τ, p) \in C (\bar{z}, τ, p) \cap N} γ_{m} (\bar{z}, τ, p) = 0

for all m = 1,2,…, where

γ_{m} (\bar{z}, τ, p)

is the m-th crossing number of

(\bar{z}, τ, p)

if

m \in J (\bar{z}, τ, p)

, or it is zero if otherwise.

Obviously, if (ii) of Lemma 5.1 is not true, then $C (\bar{z}, τ, p)$ is unbounded. Thus, if the projections of $C (\bar{z}, τ, p)$ onto z-space and onto p-space are bounded, then the projections of $C (\bar{z}, τ, p)$ onto τ-space is unbounded. Further, if we can show that the projections of $C (\bar{z}, τ, p)$ onto τ-space is away from zero, then the projections of $C (\bar{z}, τ, p)$ onto τ-space must include the interval $[τ, \infty)$ . Based on the idea, we can prove our results on the global continuation of the local Hopf bifurcation.

Lemma 5.2. If τ is bounded, then all periodic solutions to (1.6) is uniformly bounded.

Proof. Let $(x_{1} (t), x_{2} (t))$ be a nontrivial solution to (1.6) and define

x_{1} (ξ_{1}) = \min \{x_{1} (t)\}, x_{1} (η_{1}) = \max \{x_{1} (t)\}, x_{2} (ξ_{2}) = \min \{x_{2} (t)\}, x_{2} (η_{2}) = \max \{x_{2} (t)\}

with initial value

x_{1} (t) = φ_{1} (t) \geq 0, x_{1} (0) = φ_{1} (0) > 0; x_{2} (t) = φ_{2} (t) \geq 0, x_{2} (0) = φ_{2} (0) > 0,

where

t \in [- τ, 0]

. Then, it follows from (1.6) that

\{\begin{array}{c} x_{1} (t) = x_{1} (0) \exp \{\int_{0}^{t} [a - b x_{1} (t - τ) - \frac{c x_{2} (t - τ)}{m x_{2} (t - τ) + x_{1} (t - τ)}] ds\}, \\ x_{2} (t) = x_{2} (0) \exp \{\int_{0}^{t} [- d + \frac{r x_{1} (t - τ)}{m x_{2} (t - τ) + x_{1} (t - τ)}] ds\} \end{array}

Thus, $x_{1} (t) > 0, x_{2} (t) > 0,$ which implies that solutions to (1.6) cannot cross the x-axes and y-axes. Thus, the nonconstant periodic orbits must be located in the interior of the first quadrant. If (x₁(t), x₂(t)) is a solution to (1.6), then it follows from the first equation of (1.6) that

\frac{d x_{1}}{d t} < x_{1} [a - b x_{1} (t - τ)]

(5.2)

Clearly, $x_{1} (t) < x_{1} (t - τ) e^{a τ}$ for t > τ, which implies that $x_{1} (t - τ) > x_{1} (t) e^{- a τ}$ . This, together with (5.2), leads to

\frac{d x_{1}}{d t} < x_{1} [a - b e^{- a τ} x_{1} (t)],

for t > τ. Then, we obtain that

\lim_{t \to \infty} \sup x_{1} (t) \leq \frac{a e^{a τ}}{b}

Thus,

x_{1} (t) \leq \frac{a e^{a τ}}{b} + ε ≔ M

From the second equation of (1.6), we have

\frac{d x_{2}}{d t} < (r - d) x_{2},

which leads to

x_{2} (t) < e^{(r - d) τ} x_{2} (t - τ)

Then,

x_{2} (t - τ) > e^{(d - r) τ} x_{2} (t)

Therefore,

x_{2} (η_{2} - τ) > e^{(d - r) τ} x_{2} (η_{2})

Applying the second equation of (1.6), we get

- d + \frac{r x_{1} (η_{2} - τ)}{m x_{2} (η_{2} - τ) + x_{1} (η_{2} - τ)} = 0,

i.e.

\frac{r x_{1} (η_{2} - τ)}{m x_{2} (η_{2} - τ) + x_{1} (η_{2} - τ)} = d

It follows that

\frac{r M}{m e^{(d - r) τ} x_{2} (η_{2}) + M} > d,

which leads to

x_{2} (η_{2}) < \frac{(r - d) M}{d m e^{(d - r) τ}} ≔ N

Thus, the possible periodic solutions lying in the first quadrant of (1.6) must be uniformly bounded. This completes the proof of Lemma 5.2.

Lemma 5.3.If a < d, r < c, then system (1.6) has nonconstant periodic solution with period τ.

Proof. Suppose for a contradiction that system (1.6) has a nonconstant periodic solution with period τ. Then, the following ordinary differential equations

\{\begin{array}{c} \frac{d x_{1}}{d t} = x_{1} [a - b x_{1} (t) - \frac{c x_{2} (t)}{m x_{2} (t) + x_{1} (t)}], \\ \frac{d x_{2}}{d t} = x_{2} [- d + \frac{r x_{1} (t)}{m x_{2} (t) + x_{1} (t)}] \end{array}

(5.3)

has nonconstant periodic solution. System (5.3) has a boundary equilibrium E₀(0, 0) and a positive equilibrium

E_{*} (x_{1}^{*}, x_{2}^{*})

, where

x_{1}^{*} = \frac{a m^{2} + c d - c r}{a b m}, x_{2}^{*} = \frac{(r - d) (a m^{2} + c d - c r)}{a b d m^{2}}

Note that x-axis and y-axis are the invariable manifold of system (5.3) and the orbits of system (5.3) do not intersect each other. Thus, there are no solutions crossing the coordinate axes. On the other hand, consider that if system (5.3) has a periodic solution, then there must be an equilibrium in its interior, and that E₀(0, 0) is located on the coordinate axis. Thus, we can conclude that the periodic orbit of system (5.3) must lie in the first quadrant.

In the sequel, we will prove that system (5.3) has no nonconstant periodic solution in the first quadrant.

Define

D = \{(x_{1}, x_{2}) \in R^{2} | 0 \leq x_{1} \leq M, 0 \leq x_{2} \leq N\} .

It is easy to show that D is an ultimately bounded region (or absorbing and positively invariant set) of system (5.3). Let

\begin{array}{l} P (x, y) = x_{1} [a - b x_{1} (t) - \frac{c x_{2} (t)}{m x_{2} (t) + x_{1} (t)}], \\ Q (x, y) = x_{2} [- d + \frac{r x_{1} (t)}{m x_{2} (t) + x_{1} (t)}] \end{array}

Then, a direct computation show that, for (x, y) ∈ D and a < d, r < c,

\frac{\partial P (x, y)}{\partial x} + \frac{\partial Q (x, y)}{\partial y} = a - d - 2 b x_{1} + \frac{(r - c) m x_{2}^{2}}{{(m x_{2} + x_{1})}^{2}} < 0

Thus, the Bendixson–Dulac criterion, together with the fact that D is an ultimately bounded region of (5.3), implies that (5.3) has no nontrivial periodic solutions, leading to a contradiction. Thus, the proof is complete.

Theorem 5.4.Assume that a < d, r < c. Let ω₀ and τ_k (k = 0,1,2,…) be defined by (2.8) and (2.9), respectively. Then, for each τ > τ_k (k ≥ 1), system (1.6) has at least k periodic solutions.

Proof. Obviously, $(E_{*}, τ_{k}, \frac{2 π}{ω_{0}})$ is an isolated center of (1.6). Let

$C (E_{*}, τ_{k}, \frac{2 π}{ω_{0}})$ denote the connected component passing through $(E_{*}, τ_{k}, \frac{2 π}{ω_{0}})$ in Γ. It follows from Theorem 2.4 that $C (E_{*}, τ_{k}, \frac{2 π}{ω_{0}})$ is nonempty. It is only necessary to prove that the projection of $C (E_{*}, τ_{k}, \frac{2 π}{ω_{0}})$ onto τ-space is $[\bar{τ}, \infty)$ for each j > 0, where $\bar{τ} \leq τ_{k}$ .

By Lemma 2.2 and 2.3, there exists ɛ > 0, δ > 0 and a smooth curve λ: (τ_k – δ, τ_k + δ) → C such that $det Δ (λ (τ)) = 0, | λ (τ) - i ω_{0} | < ε$ for all $τ \in [τ_{k} - δ, τ_{k} + δ]$ , and $λ (τ_{k}) = i ω_{0}, \frac{d}{d τ} Re (λ (τ)) |_{τ = τ_{k}} > 0$

Let $Ω_{ε, \frac{2 π}{ω_{0}}} = \{(μ, p) : 0 < μ < ε, | p - \frac{2 π}{ω_{0}} | < ε\}$ . It is easy to verify that on $[τ_{k} - δ, τ_{k} + δ] \times \partial Ω_{ε, \frac{2 π}{ω_{0}}}$ , $det Δ_{(E_{*}, τ_{k}, p)} (μ + \frac{2 π}{p} i) = 0$ if and only if μ = 0, τ = τ_k and $p = \frac{2 π}{ω_{0}}$ . This verifies assumption (A₄) of Wu [29]. Moreover, if we put

H^{\pm} (E_{*}, τ_{k}, \frac{2 π}{ω_{0}}) = det (Δ_{(E_{*}, τ_{k}, p)} (μ + i \frac{2 π}{p})),

then, we have the crossing number of isolated centers

(E_{*}, τ_{k}, \frac{2 π}{ω_{0}})

as follows:

γ_{1} (E_{*}, τ_{k}, \frac{2 π}{ω_{0}}) = {deg}_{B} (H^{-} (E_{*}, τ_{k}, \frac{2 π}{ω_{0}}), Ω_{ε, \frac{2 π}{ω_{0}}}) - {deg}_{B} (H^{+} (E_{*}, τ_{k}, \frac{2 π}{ω_{0}}), Ω_{ε, \frac{2 π}{ω_{0}}}) = - 1

By Theorem 3.3 of Wu [29], we conclude that the connected component $C (E_{*}, τ_{k}, \frac{2 π}{ω_{0}})$ through $(E_{*}, τ_{k}, \frac{2 π}{ω_{0}})$ in Σ is nonempty, and

\sum_{(\bar{z}, τ, p) \in C (\bar{z}, τ, p)} γ (\bar{z}, τ, p) < 0

Thus, $C (E_{*}, τ_{k}, \frac{2 π}{ω_{0}})$ is unbounded. From (2.9), one can show that, for k ≥ 1,

τ_{k} = \frac{1}{ω_{0}} [\arcsin \frac{(m_{1} + n_{2}) \sqrt{m_{1} n_{2} - m_{2} n_{1}}}{2 (m_{1} n_{2} - m_{2} n_{1})} + 2 k π], k = 1, 2, \dots

Hence,

\frac{2 π}{ω_{0}} < τ_{k}

(5.4)

Now, we are in position to prove that the projection of $C (E_{*}, τ_{k}, \frac{2 π}{ω_{0}})$ onto τ-space is [τ, ∞), where τ ≤ τ_k. Clearly, It follows from the proof of Lemma 5.3 that system (1.6) with τ = 0 has no nontrivial periodic solution. Hence, the projection of $C (E_{*}, τ_{k}, \frac{2 π}{ω_{0}})$ onto τ-space is away from zero.

For the sake of contradiction, we suppose that the projection of $C (E_{*}, τ_{k}, \frac{2 π}{ω_{0}})$ onto τ-space is bounded. This means that the projection of $C (E_{*}, τ_{k}, \frac{2 π}{ω_{0}})$ onto τ-space is included in an interval (0, τ^*). From (5.4) and applying Lemma 5.3, we get 0 < p < τ* for $(z_{t}, τ, p) \in C (E_{*}, τ_{k}, \frac{2 π}{ω_{0}})$ . This implies that the projection of $C (E_{*}, τ_{k}, \frac{2 π}{ω_{0}})$ onto p-space is also bounded. Thus, combining this with Lemma 5.2, we can get that the connected component $C (E_{*}, τ_{k}, \frac{2 π}{ω_{0}})$ is bounded. This contradiction completes the proof.

6 Conclusions and biological explanations

In this paper, we have investigated the local stability of the positive equilibrium $E_{*} (x_{1}^{*}, x_{2}^{*})$ and the local Hopf bifurcation in a delayed predator–prey model with Hassell–Varley-type functional response. We have showed that if conditions (H₁)–(H₃) hold, the positive equilibrium $E_{*} (x_{1}^{*}, x_{2}^{*})$ of system (1.6) is asymptotically stable for all τ ∈ [0, τ₀) and unstable for τ > τ₀. This shows that, in this case, the population of preys and predators will tend to stabilization and still keep stable whenever the delay parameter lies in the range τ ∈ [0, τ₀) and unstable when τ > τ₀. We have also showed that, if conditions (H₁)–(H₃) hold, as the delay τ increases, the equilibrium loses its stability and a sequence of Hopf bifurcations occur at $E_{*} (x_{1}^{*}, x_{2}^{*})$ , i.e., a family of periodic orbits bifurcate from the positive equilibrium $E_{*} (x_{1}^{*}, x_{2}^{*})$ . This means that the population of preys and predators may coexist and keep in an oscillatory mode. Moreover, the direction of Hopf bifurcation and the stability of the bifurcating periodic orbits are discussed by applying the normal form theory and the center manifold theorem. Numerical simulations supporting our theoretical results are also included. Further, sufficient conditions ensuring the existence of global Hopf bifurcation are given, i.e., if a < d, r < c, then system (1.6) has at least k periodic solutions for τ > τ_k(k ≥ 1). It is shown that the population of preys and predators still keep in an oscillatory mode near the positive equilibrium $E_{*} (x_{1}^{*}, x_{2}^{*})$ for τ > τ_k(k ≥ 1).

References

[1] V. Volterra Variazioni e fluttuazioni del numero di individui in specie animali conviventi, Mem. Accd. Lincei., Volume 2 (1926), pp. 31-113

[2] A. Lotka Elements of Physical Biology, Williams and Wilkins, Baltimore, MD, USA, 1925

[3] R. Arditi; L.R. Ginzburg Coupling in predator–prey dynamics: ratio-dependence, J. Theoret. Biol., Volume 139 (1989), pp. 311-326

[4] H.K. Baek Qualitative analysis of Beddington-DeAngelis type impulsive predator–prey models, Nonlinear Anal.: Real World Appl., Volume 11 (2010), pp. 1312-2132

[5] H. Baek; Y. Lim Dynamics of an impulsively controlled Michaelis–Menten-type predator–prey system, Commun. Nonlinear Sci. Numer. Simul., Volume 16 (2011), pp. 2041-2053

[6] Y.L. Song; S.L. Yuan; J.M. Zhang Bifurcation analysis in the delayed Leslie–Gower predator–prey system, Appl. Math. Model., Volume 33 (2009), pp. 4049-4061

[7] C.J. Xu; X.H. Tang; M.X. Liao Stability and bifurcation analysis of a delayed predator–prey model of prey dispersal in two-patch environments, Appl. Math. Comput., Volume 216 (2010), pp. 2920-2936

[8] C.J. Xu; M.X. Liao; X.F. He Stability and Hopf bifurcation analysis for a Lokta–Volterra predator–prey model with two delays, Int. J. Appl. Math. Comput. Sci., Volume 21 (2011), pp. 97-107

[9] C.J. Xu; X.H. Tang; M.X. Liao; X.F. He Bifurcation analysis in a delayed Lokta–Volterra predator–prey model with two delays, Nonlinear Dynam., Volume 66 (2011), pp. 168-183

[10] J.Z. Zhang; J. Zhen; J. Yan; G.Q. Sun Stability and Hopf bifurcation in a delayed competition system, Nonlinear Anal.: Real World Appl., Volume 70 (2009), pp. 658-670

[11] T. Zhao; Y. Kuang; H.L. Simith Global existence of periodic solutions in a class of delayed Gause-type predator–prey system, Nonlinear Anal.: Real World Appl., Volume 28 (1997), pp. 1373-1394

[12] B. Rai; W. Krawcewicz Hopf bifurcation in symmetric configuration of predator–prey-mutualist systems, Nonlinear Anal.: Real World Appl., Volume 71 (2009), pp. 4279-4296

[13] N. Bairagi; D. Jana On the stability and Hopf bifurcation of a delay-induced predator–prey system with habitat complexity, Appl. Math. Modell., Volume 35 (2011), pp. 3255-3267

[14] M. Baurmann; T. Gross; U. Feudel Instabilities in spatially extended predator–prey systems: spatio-temporal patterns in the neighborhood of Turing-Hopfbifurcations, J. Theor. Biol., Volume 245 (2007), pp. 220-229

[15] S. Kumar; S.K. Srivastava; P. Chingakham Hopf bifurcation and stability analysis in a harvested one-predator-two-prey model, Appl. Math. Comput., Volume 129 (2002), pp. 107-118

[16] T.K. Kar; A. Batabyal Stability and bifurcation of a prey–predator model with time delay, C.R. Biol., Volume 332 (2009), pp. 642-651

[17] C. Çelik Hopf bifurcation of a ratio-dependent predator–prey system with time delay, Chaos, Solitons Fractals, Volume 42 (2009), pp. 1474-1484

[18] C. Çelik The stability and Hopf bifurcation for a predator–prey system with time delay, Chaos, Solitons Fractals, Volume 37 (2008), pp. 87-99

[19] K. Chakraborty; M. Chakraborty; T.K. Kar Optimal control of harvest and bifurcation of a prey–predator model with stage structure, Appl. Math. Comput., Volume 217 (2011), pp. 8778-8792

[20] I. Kusbeyzi; O.O. Aybar; A. Hacinliyan Stability and bifurcation in two species predator–prey models, Nonlinear Anal.: Real World Appl., Volume 12 (2011), pp. 377-387

[21] P. Auger; B.W. Kooi; R.B. de la Parra; J.C. Poggiale Bifurcation analysis of a predator–prey model with predators using hawk and dove tactics, J. Theoret. Biol., Volume 238 (2006), pp. 597-607

[22] S. Kovács; K. Kiss; M. Farkas Qualitative behaviour of a ratio-dependent predator–prey system, Nonlinear Anal.: Real World Appl., Volume 10 (2009), pp. 1627-1642

[23] X.P. Yan; C.H. Zhang Hopf bifurcation in a delayed Lokta–Volterra predator–prey system, Nonlinear Anal.: Real World Appl., Volume 9 (2008), pp. 114-127

[24] S.L. Yuan; Y.L. Song Stability and Hopf bifurcations in a delayed Lesile-Gower predator–prey system, J. Math. Anal. Appl., Volume 355 (2009), pp. 82-100

[25] H. Freedman Deterministic Mathematical Models in Population Ecology, Marcel Dekker, New York, 1980

[26] M. Hassell; G. Varley New inductive population model for insect parasites and its bearing on biological control, Nature, Volume 223 (1969), pp. 1133-1136

[27] K. Wang Periodic solutions to a delayed predator–prey model with Hassell–Varley-type functional response, Nonlinear Anal.: Real World Appl., Volume 12 (2011), pp. 137-145

[28] T. Faria; L.T. Magalhaes Normal forms for retared functional differential equations with parameters and applications to Hopf bifurcations, J. Differ. Equations, Volume 122 (1995), pp. 181-200

[29] J. Wu Symmertic functional differential equations and neural networks with memory, Trans. Am. Math. Soc., Volume 350 (1998), pp. 4799-4838

[30] S. Ruan; J. Wei On the zero of some transcendential functions with applications to stability of delay differential equations with two delays, Dyn. Contin. Discrete Impuls. Syst. Ser. A, Volume 10 (2003), pp. 863-874

[31] X. Dieuonnné Founcations of Modern Analysis, Academic Press, New York, 1960

[32] J.K. Hale Theory of Functional Differential Equations, Springer-Verlag, New York, 1977

[33] S.N. Chow; J.K. Hale Methods of Bifurcation Theory, Springer, New York, 1982

^☆ This work is supported by National Natural Science Foundation of China (No. 11261010 and No. 60902044), the Doctoral Foundation of Guizhou College of Finance and Economics (2010), the Science and Technology Program of Hunan Province (No. 2010FJ6021) and the soft Science and Technology Program of Guizhou Province (No. 2011LKC2030).

Comments - Policy

1 Introduction

H1

2 Stability of the equilibrium and existence of the local Hopf bifurcation

H2

H3

3 Direction and stability of the Hopf bifurcation

4 Numerical examples

5 Global continuation of local Hopf bifurcations

6 Conclusions and biological explanations

H₁

H₂

H₃