Hopf bifurcation in a delayed predator-prey system with asymmetric functional response and additional food

1.   Introduction

In the real world, we often need to control or eliminate the target prey in a prey-predator system. Biological control and chemical control are the two essential methods. The main means of chemical control is to spray pesticides on pests at different fixed times (see e.g. ^[1]). However, chemical control is said to be environmentally detrimental. There are various forms of biological control, which is the focus of current research, such as immunocontraceptive technology, introducing other predators, offering the predator additional food and so on. Saunders et al. ^[2] used the approaches of immunocontraception and daughterless genes to control the growth of the target pest species. Ghosh et al. ^[3] investigated the effectiveness of periodic impulsive releases of natural enemies into a two-patch environment.

Numerous studies had shown that the additional food (to predator) can be helpful to increase predators and the prey can be controlled. Tena et al. ^[4] showed that the melinus females were better able to forage and oviposit by providing sugar to them. Srinivasu et al. ^[5] studied a standard predator-prey model with additional food and presented the evidence that how to eliminate target pests. Sahoo et al. ^[6] investigated that the chaotic population dynamics can be controlled to obtain regular population dynamics only by supplying additional food to top predator. Research on the models with additional food, one can refer to the literature ^{[7,8,9,10,11]}.

Recently, Basheer et al. ^[12] believed that the additional food can increase the carrying capacity and birth rate of the predators, so they studied the following Holling-Tanner model with additional food:

where $u(t)$ and $v(t)$ denote the prey and predator density. $\alpha$, $\beta\in \mathbb{R}_+$, and $\frac{1} {\alpha}$ and $\beta$ represent the quality and quantity of the additional food, respectively. The biological significance of other parameters refer to ^[12]. Their research indicated that a conditional stable prey-extinction equilibrium could be obtained in system (1.1), but it was nonexistent in the absence of additional food.

On the other hand, if the population maturation time is considered, the corresponding model should be delayed. Therefore, lots of delayed predator-prey systems were studied (see, e.g. ^{[13,14,15,16,17,18,19,20,21,22,23,24]}). It is interesting that Jiang et al. ^[25] applied a delay differential equation (DDE) to two-enterprise interaction mechanism. The dynamical properties of DDE are far more complex than those of ordinary differential equation (ODE). A time delay can cause an equilibrium undergoes from stable state to unstable, which leads to the complicated dynamics of DDE. Arising periodic solution through the Hopf bifurcation is one of the hot topics (see, e.g. ^{[18,19,20,21,22,23,24,26]}). However, these bifurcating periodic solutions from Hopf bifurcation are generally local. Whether these local periodic solutions exist globally is an interesting subject. Erbe et al. ^[27] established the global Hopf bifurcation theorem and Wu ^[28] applied it to a neural networks with memory. Thereafter, some researchers employed the theorem in ^[28] to study the global existence of periodic solutions for DDE (see, e.g. ^{[29,30,31,32,33,34]}).

Inspiration based on model (1.1), we consider the maturation time of predator and standard Holling type II functional response in the model, then we improve model (1.1) as follows:

where $\tau$ represent the maturation time of predator. The initial conditions are chosen as:

where $(\varphi_1(\theta), \varphi_2(\theta))\in C\{[-\tau, 0], \mathbb{R}_+^2\}, \mathbb{R}_+^2 = \{u, v:u\geq 0, v\geq 0\}$.

The main purpose of our work are concluded as follows:

(I) We consider the feedback time delay in the maturation time of predator, which is more general than the works in ^[12]. We investigate the effects of feedback time on the stability of the equilibria and the conditions on occurring Hopf bifurcations.

(II) The delayed feedback is designed for studying bifurcating periodic solutions. We analyze the direction and stability of bifurcating periodic solutions on the center manifold.

(III) The global existence of bifurcating periodic solutions are studied mathematically.

The rest of this paper is organized as follows. We first discuss some properties of system (1.2) to prepare for the next section. In Section 3, we study the local stability of each feasible equilibrium of system (1.2) with the effect of the time delay $\tau$. The formulas determining the direction and stability of bifurcating periodic solutions are obtained via the theory in Hassard et al. ^[35] in Section 4. In Section 5, the global existence of periodic solutions under the second critical value is proved by using the theory in Wu ^[28]. Finally, some examples are utilized to demonstrate the validity of the previous results.

2.   Preliminaries

Clearly, system (1.2) can be calculated by

According to positive initial values of system (1.2), it is not difficult to obtain the following lemma.

Lemma 2.1 Any of the solutions of system (1.2) are positive for $t\geq0$ with positive initial values.

Theorem 2.2 System (1.2) is ultimately bounded when $\tau$ is bounded.

Proof. From the first equation of system (1.2) and by using the comparison theorem, it is easy to obtain

Namely, there exists a time $T_1$ such that $u(t)\leq 1+\varepsilon$ for arbitrary $\varepsilon > 0$ and $t > T_1$. By the second equation of (1.2), we have

Integrating both sides on the interval $[t-\tau, t]$, it produces

which implies

Meanwhile, it follows from the second equation of (1.2) that

Using the comparison theorem, we obtain

Therefore, Theorem 2.2 is confirmed.

Theorem 2.3 In the absence of delay, if $s < \delta$, then there is no closed loop in the first quadrant of system (1.2).

Proof. Let $B(u, v) = \frac{1}{uv}$, $f_1 = u(1-u)-\frac{suv}{m+\alpha\beta+u}$, $f_2 = \delta v\Big(n+\frac{\beta-v}{\alpha\beta+u)}\Big)$, then

By using the Bendixson-Dulac criterion, the proof is completed.

3.   Equilibrium, local stability and Hopf bifurcation

In the paper, we use the following representations for the sake of simplicity:

where $(\bar{u}, \bar{v})$ stands for the positive interior equilibrium and is defined in this section.

In order to obtain the equilibria, we discuss the algebraic equations:

Let $(\bar{u}, \bar{v})$ stands for the interior equilibrium, where $\bar{u}$ is the positive root of the equation $u^2+Au+B = 0$. Thus, we have $\bar{u}_\pm = \frac{-A\pm\sqrt{A^2-4B}}{2}$ and $\bar{v} = n(\alpha\beta+\bar{u})+\beta$. The equilibria are as follows:

(i) Trivial equilibrium $E_0 = (0, 0)$.

(ii) Predator-extinction equilibrium $E_1 = (1, 0)$.

(iii) Prey-extinction equilibrium $E_2 = (0, n\alpha\beta+\beta)$.

(iv) A unique coexisting equilibrium $E_+ = (\bar{u}_+, \bar{v}_+)$ when $B < 0$; two coexisting equilibria $E_\pm = (\bar{u}_\pm, \bar{v}_\pm)$ when $B > 0, \ A < 0$ and $A^2-4B > 0$.

Let $E = (u^*, v^*)$ be arbitrary equilibrium. We use linearization technique to analyze the local stability of system (1.2). The Jacobian matrix of system (1.2) at $E = (u^*, v^*)$ is given by

It is easy to confirm that $E_0 = (0, 0)$ is an unstable node and $E_1 = (1, 0)$ is a saddle.

At Prey-extinction equilibrium $E_2 = (0, n\alpha\beta+\beta)$, the corresponding Jacobian matrix is

and the characteristic equation becomes

where $P_1 = 1-\frac{s(\alpha\beta n+\beta)}{\alpha\beta}$, $P_2 = \delta(n+\frac{1}{\alpha})$.

Case 1. $\tau = 0$.

It follows from (3.2) that $E_2 = (0, n\alpha\beta+\beta)$ is locally asymptotically stable when $P_1 < 0$ (equivalent to $B > 0$) and unstable when $P_1 > 0$.

Case 2. $\tau > 0$.

Let $\lambda = i\omega^* (\omega^* > 0)$ be a root of the equation $\lambda+P_2e^{-\lambda\tau} = 0$, then $i\omega+P_2(\cos\omega\tau-i\sin\omega\tau) = 0$. By a direct calculation, we get $\omega_0^* = P_2$ and $\tau_k^* = \frac{1}{P_2}(2k\pi+\frac{\pi}{2}), k = 0, 1, 2, \cdots$. Differentiating the both sides of $\lambda+P_2e^{-\lambda\tau} = 0$ with respect to $\tau$, we have

that is

which implies that

Lemma 3.1 If $B > 0$, then all roots of the characteristic equation (3.2) have negative real part when $0\leq\tau < \frac{\pi}{2P_2}$ and at least one positive real part when $\tau > \frac{\pi}{2P_2}$.

Therefore, we have the following conclusions for the boundary equilibria.

Theorem 3.2 (i) The trivial equilibrium $E_0 = (0, 0)$ and predator-extinction equilibrium $E_1 = (1, 0)$ are always unstable for all $\tau\geq0$.

(ii) When $B > 0$, the prey-extinction equilibrium $E_2 = (0, n\alpha\beta+\beta)$ is asymptotically stable for all $0\leq\tau < \frac{\pi}{2P_2}$ and unstable for all $\tau > \frac{\pi}{2P_2}$. The system (1.2) undergoes a Hopf bifurcation at $E_2$ for $\tau = \frac{\pi}{2P_2}$.

Now we investigate the stability of the coexisting equilibrium $E = (\bar{u}, \bar{v})$. The corresponding Jacobian matrix is

and the corresponding characteristic equation becomes

Case 1. $\tau = 0$.

Equation (3.3) becomes

We calculate $F$ as follows:

which implies that $F < 0$ when $\bar{u} = \bar{u}_-$ and $F > 0$ when $\bar{u} = \bar{u}_+$. Thus, the following conclusions are obvious.

Lemma 3.3 If the coexisting equilibria exist and $\tau = 0$, then $E_+ = (\bar{u}_+, \bar{v}_+)$ is locally asymptotically stable when $G+D > 0$ and $E_- = (\bar{u}_-, \bar{v}_-)$ is always unstable.

Case 2. $\tau > 0$.

Let $\lambda = i\omega (\omega > 0)$ be a root of the equation (3.3), then

We obtain

that is

From (3.6), we have

Let $z = \omega^2$, (3.7) turns to

Clearly, (3.8) has a unique positive root $z_0 = \frac{-(G^2-D^2)+\sqrt{(G^2-D^2)^2+4D^2F^2}}{2}$. So (3.8) has a pair of purely imaginary roots $\pm i\omega_0$($\omega_0 = \sqrt{z_0}$). From the first equation of (3.6), we have

Differentiating both sides of (3.3) with respect to $\tau$, we have

namely,

Then

which implies that

Lemma 3.4 If condition $G+D > 0$ holds, then all roots of the characteristic equation (3.3) at $E_+ = (\bar{u}_+, \bar{v}_+)$ have negative real part when $0 < \tau < \tau_0$ and at least one positive real part when $\tau > \tau_0$.

From lemma 3.3 and 3.4, we confirm the following conclusions about the interior equilibrium.

Theorem 3.5 (i) If $E_- = (\bar{u}_-, \bar{v}_-)$ exists, then it is unstable for all $\tau\geq0$.

(ii) If $E_+ = (\bar{u}_+, \bar{v}_+)$ exists and the condition $G+D > 0$ holds, then $E_+$ is locally asymptotically stable for all $0\leq\tau < \tau_0$ and unstable for all $\tau > \tau_0$. The system (1.2) undergoes a Hopf bifurcation at $E_+$ for $\tau = \tau_0$.

Remark 3.6 If the prey is pest, we just need to control the quantity of additional food satisfy the following conditions:

then additional food can induce pest eradication. Meanwhile, the density of predators eventually goes to $n\alpha\beta+\beta$ when the maturation time of predator species is less than $\frac{\pi}{2P_2}$.

Remark 3.7 If the conditions $B > 0$, $A < 0$, $A^2-4B > 0$ and $G+D > 0$ hold simultaneously, the system (1.2) is bistable when $\tau < \text{min}\{\tau_0^*, \tau_0\}$.

4.   Direction and stability of the Hopf bifurcation

We know from the literature ^[35] that the properties of Hopf bifurcation are determined by the following three quantities, namely

where $c_1(0) = \frac{i}{2\tilde{\omega}\tilde{\tau}}\Big[g_{11}g_{20}-2|g_{11}|^2-\frac{\mid g_{02}\mid ^2}{3}\Big]+\frac{g_{21}}{2}$, $\tilde{\tau}$ is the critical value. $\mu_2$, $\beta_2$ and $T_2$ determine direction, stability and the period of the bifurcating periodic solutions, respectively. So we need to figure out the value of $g_{ij}$ in $c_1(0)$. We will use the normal form theory and the center manifold theorem to obtain the expression of $g_{ij}$ in this section.

It has known that system (1.2) undergoes Hopf bifurcation at coexisting equilibrium $E_+$ when $\tau = \tau_0$. We denote the critical values $\tau_k$ and $E_+ = (\bar{u}_+, \bar{v}_+)$ as $\tilde{\tau}$ and $E^* = (u^*, v^*)$, respectively. Let $x(t) = u(\tau t)-u^\ast$ and $y(t) = v(\tau t)-v^\ast$, using Taylor expansion, (1.2) can be rewritten as

where

$a_1 = 1-2u^\ast-\frac{s(m+\alpha\beta)v^*}{(m+\alpha\beta+u^\ast)^2}$, $a_2 = -\frac{su^\ast}{m+\alpha\beta+u^\ast}$, $a_3 = -\frac{\delta\nu^\ast(\beta-v^*)}{(\alpha\beta+u^\ast)^2} = \frac{\delta nv^*}{\alpha\beta+u^\ast}$,

$a_4 = -\frac{\delta v^*}{\alpha\beta+u^\ast}$, $a_5 = -2+\frac{2s(m+\alpha\beta)v^*}{(m+\alpha\beta+u^\ast)^3}$, $a_6 = -\frac{2s(m+\alpha\beta)}{(m+\alpha\beta+u^\ast)^2}$,

$a_7 = \frac{2\delta\nu^\ast(\beta-v^*)}{(\alpha\beta+u^\ast)^3}$, $a_8 = -\frac{2\delta(\beta-v^*)}{(\alpha\beta+u^\ast)^2}$, $a_9 = \frac{2\delta v^*}{(\alpha\beta+u^\ast)^2}$, $a_{10} = -\frac{2\delta}{\alpha\beta+u^\ast}$.

Let $\tau = \bar{\tau}+h$, then $h = 0$ is a Hopf bifurcation value of the system (1.2). Choose the phase space $C = C([-1, 0], \mathbb{R}^2)$, $\phi(\theta) = (\phi_1(\theta), \phi_2(\theta))^T\in C, \theta\in[-1, 0]$, define $L(h)$

By the Riesz representation theorem, we choose the bounded variation function

such that

where $\delta(\theta)$ is delta function.

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

and

where $T_0(\theta) = \left\{ \begin{aligned} \; I, & \ \; \; \theta = 0, \\ \; 0, & \ \; \; \theta\in[-1, 0). \end{aligned} \right.$

Then (4.2) is written as

where $u_t = u(t+\theta)$ and $u_t = (x_t, y_t)^\text{T}$.

For $\phi\in \mathbb{C}^1([-1, 0], \mathbb{C}^2)$ and $\psi\in \mathbb{C}^1([0, 1], (\mathbb{C}^2)^*)$, define a adjoint operator of $A(0)$

and the bilinear form

From the previous discussion, we know $\pm i\tilde{\omega}\tilde{\tau}$ are the eigenvalues of $A(0)$ and $A^*$. Let $q(\theta) = (1, \alpha_1)^T e^{i\tilde{\omega}\tilde{\tau}\theta}$ be the eigenvector of $A(0)$ corresponding to $i\tilde{\omega}\tilde{\tau}$ and $q^*(s) = M(1, \alpha_2) e^{i\tilde{\omega}\tilde{\tau}s}$ be the eigenvector of $A^*$ corresponding to $-i\tilde{\omega}\tilde{\tau}$, and they satisfy the conditions $\langle q^*, q\rangle = 1$ and $\langle q^*, q\bar{}\rangle = 0$. Therefore, we have

and

then $\alpha_1 = \frac{i\tilde{\omega}-a_1}{a_2}$ and $\alpha_2 = \frac{-i\tilde{\omega}-a_1}{a_3e^{i\tilde{\omega}\tilde{\tau}}}$.

By the bilinear form, we have

Thus $M = \frac{1}{(1+\bar{\alpha}_1\alpha_2)+\tilde{\tau}\alpha_2(a_3+a_4\bar{\alpha}_1)e^{i\tilde{\omega}\tilde{\tau}}}$.

We need the coordinates to describe the center manifold $C_0$ near $h = 0$. Let $z$ and $\bar{z}$ be local coordinates for $C_0$ in the directions of $q^*$ and $\bar{q}^*$. Assume that $u_t$ is a solution of (4.3) at $h = 0$, define

On $C_0$, we have $W(t, \theta) = W(z(t), \bar{z}(t), \theta)$, where

The manifold of (4.2) on the center manifold is determined by the following equation

which is abbreviated to

where the power series form of $g(z, \bar{z})$ is

and

By a direct calculation, we have

Substituting (4.8) into $F_0$ and comparing with (4.7), we obtain

$g_{20} = 2\tilde{\tau}\bar{M}(a_5+a_6\alpha_1+\bar{\alpha}_2(a_7e^{-2i\tilde{\omega}\tilde{\tau}}+ a_8\alpha_1e^{-i\tilde{\omega}\tilde{\tau}}+a_9\alpha_1e^{-2i\tilde{\omega}\tilde{\tau}} +a_{10}\alpha_1^2e^{-i\tilde{\omega}\tilde{\tau}}))$,

$g_{02} = 2\tilde{\tau}\bar{M}(a_5+a_6\bar{\alpha}_1+\bar{\alpha}_2(a_7e^{2i\tilde{\omega}\tilde{\tau}}+ a_8\bar{\alpha}_1e^{i\tilde{\omega}\tilde{\tau}}+a_9\bar{\alpha}_1e^{2i\tilde{\omega}\tilde{\tau}} +a_{10}\bar{\alpha}_1^2e^{i\tilde{\omega}\tilde{\tau}}))$,

$g_{11} = 2\tilde{\tau}\bar{M}(a_5+a_6Re\{\alpha_1\}+\bar{\alpha}_2(a_7+a_8Re\{\alpha_1e^{i\tilde{\omega}\tilde{\tau}}\} +a_9Re\{\alpha_1\}+a_{10}Re\{\alpha_1\bar{\alpha}_1e^{i\tilde{\omega}\tilde{\tau}}\}))$,

$g_{21} = 2\tilde{\tau}\bar{M}(a_5k_5+a_6k_6+\bar{\alpha}_2(a_7k_7+a_8k_8+a_9k_9+a_{10}k_{10}))$,

where

$k_5 = W_{20}^{(1)}(0)+2W_{11}^{(1)}(0)$,

$k_6 = \frac{1}{2}\bar{\alpha}_1W_{20}^{(1)}(0)+\alpha_1W_{11}^{(1)}(0)+\frac{1}{2}W_{20}^{(2)}(0)+W_{11}^{(2)}(0)$,

$k_7 = e^{i\tilde{\omega}\tilde{\tau}}W_{20}^{(1)}(-1)+2e^{-i\tilde{\omega}\tilde{\tau}}W_{11}^{(1)}(-1)$,

$k_8 = \frac{1}{2}\bar{\alpha}_1W_{20}^{(1)}(-1)+\alpha_1W_{11}^{(1)}(-1)+\frac{1}{2}e^{i\tilde{\omega}\tilde{\tau}}W_{20}^{(2)}(0)+e^{-i\tilde{\omega}\tilde{\tau}}W_{11}^{(2)}(0)$,

$k_9 = \frac{1}{2}e^{i\tilde{\omega}\tilde{\tau}}W_{20}^{(2)}(-1)+e^{-i\tilde{\omega}\tilde{\tau}}W_{11}^{(2)}(-1) +\frac{1}{2}\bar{\alpha}_1e^{i\tilde{\omega}\tilde{\tau}}W_{20}^{(1)}(-1)+\alpha_1e^{-i\tilde{\omega}\tilde{\tau}}W_{11}^{(1)}(-1)$,

$k_{10} = \frac{1}{2}\bar{\alpha}_1W_{20}^{(2)}(-1)+\alpha_1W_{11}^{(2)}(-1) +\frac{1}{2}\bar{\alpha}_1e^{i\tilde{\omega}\tilde{\tau}}W_{20}^{(2)}(0)+\alpha_1e^{-i\tilde{\omega}\tilde{\tau}}W_{11}^{(2)}(0)$.

In order to obtain the normal form of (4.6) confined to the center manifold, we need compute $W_{20}(\theta)$ and $W_{11}(\theta)$.

Using $\dot{W} = \dot{u}_t-\dot{z}q-\dot{\bar{z}}\bar{q}$ which combines with (4.3) and (4.6), we obtain

On the other hand, from (4.5) and (4.6), we have

We substitute (4.6) into (4.9) and compare the coefficients of $\frac{z^2}{2}$ and $z\bar{z}$ with (4.10), respectively. It gives

and

According to (4.11) and (4.12), by a direct calculation for $\theta\in[-1, 0)$, we obatin

and

where $E_1$ and $E_2$ hold the following equations

From the above discussions, we have the following conclusions on the center manifold.

Theorem 4.1 (i) If $\mu_2 > 0$ ($\mu_2 < 0$), then the Hopf bifurcation is supercritical (subcritical).

(ii) If $\beta_2 < 0$ ($\beta_2 > 0$), the bifurcating periodic solution is stable (unstable).

(iii) If $T_2 > 0$ $(T_2 < 0)$, the period increases (decreases).

5.   Global existence of periodic solutions

Next, we investigate the global continuation of bifurcating periodic solutions from the positive equilibrium $(E_+, \tau_k)$. We employ the global Hopf bifurcation theorem and follow closely the notations in Wu ^[28].

Denote the conditions for the existence of $E_+$ as

(H) Either $B < 0$ or $B > 0, A < 0$ and $A^2-4B > 0$.

Let $\mathbb{R}_+ = \{(u, v)\in \mathbb{R}^2, u > 0, v > 0\}$, $X = C([-\tau, 0], \mathbb{R}^2)$, $z_t = (u_t, v_t)$, the system (1.2) is rewritten as

where $z_t(\theta) = z(t+\theta)\in X$ and $(\tau, p)\in \mathbb{R}_+\times \mathbb{R}_+$. Clearly, the mapping $F:X\times \mathbb{R}_+\times \mathbb{R}_+\longrightarrow \mathbb{R}^2$ is completely continuous. If we take $\mathbb{R}^2$ for the subspace of constant functions of $X$, we obtain $\hat{F}\mid_{\mathbb{R}^2\times \mathbb{R}_+\times \mathbb{R}_+}:\mathbb{R}^2\times \mathbb{R}_+\times \mathbb{R}_+\longrightarrow \mathbb{R}^2$. It is easy to know $(E_0, \tau, p)$, $(E_1, \tau, p)$, $(E_2, \tau, p)$, $(E_-, \tau, p)$ and $(E_+, \tau, p)$ are all stationary solutions of (5.1). Now, we verify that $(E_+, \tau, p)$ holds the conditions ($\text{A}_1$), ($\text{A}_2$), ($\text{A}_3$) and ($\text{A}_4$) in ^[28].

From $(1.2)$, We know easily that $\hat{F}\in C^2(\mathbb{R}_+^2\times\mathbb{R}_+\times\mathbb{R}_+, \mathbb{R}_+^2)$ and $F(\phi, \tau, p)$ is differential with respect to $\phi$. That is to say, the conditions ($\text{A}_1$) and ($\text{A}_3$) are satisfied. We also obtain

Directly calculate, we have

Then $D_z\hat{F}(E_+, \tau, p)$ is a homeomorphism on $\mathbb{R}^2$ at $E_+$, which satisfies the condition ($\text{A}_2$).

The characteristic matrix of (5.1) at $(E_+, \tau, p)$ is taken as:

A stationary solution $(E_+, \tau, p)(\lambda)$ of (5.1) is called a center if $\text{det} \Delta_{(E_+, \tau, p)}(\lambda) = 0$ has purely imaginary characteristic roots of the form $\text{i}m\frac{2\pi}{p_0}$ for some positive integer $m$. It follows from (5.2) that

which is the same as $(3.3)$. Taking $p_0 = \frac{2\pi}{\omega_0}$, we know $\text{i}\frac{2\pi}{p_0}$ is a root of $(5.3)$, namely, $\text{i}\frac{2\pi}{p_0}$ is an eigenvalue of $(E_+, \tau_k, \frac{2\pi}{\omega_0})$. Thus, $(E_+, \tau_k, \frac{2\pi}{\omega_0})$ is a center, where $\tau_k$ and $\omega_0$ are defined in Section 3. Furthermore, the center $(E_+, \tau_k, \frac{2\pi}{\omega_0})$ is an isolated center, because it satisfies the following two conclusions:

(i) $J(E_+, \tau_k, \frac{2\pi}{\omega_0}) = 1$, where $J(E_+, \tau_k, \frac{2\pi}{\omega_0})$ is a positive integer set with respect to $m$ such that $\text{i}m\frac{2\pi}{p_0}$ are eigenvalues of $(E_+, \tau_k, \frac{2\pi}{\omega_0})$.

(ii) For arbitrary $k\geq0$, there exist $\varepsilon_k > 0, \delta_k > 0$ and a smooth function $\lambda:(\tau_k-\sigma_k, \tau_k+\sigma_k\rightarrow C)$(C is complex field) such that $\text{det}\Delta_{(E_+, \tau_k, \frac{2\pi}{\omega_0})}(\lambda(\tau)) = 0$, $|\lambda(\tau)-\text{i}\omega_0| < \varepsilon_k$ for arbitrary $\tau\in(\tau_k-\sigma_k, \tau_k+\sigma_k)$ and $\lambda(\tau_k) = \text{i}\omega_0$, $\frac{\text{d}\lambda(\tau)}{\text{d}\tau}|_{\tau = \tau_k} > 0$. That is, $(E_+, \tau_k, \frac{2\pi}{\omega_0})$ is the only center in certain neighborhood of $(E_+, \tau_k, \frac{2\pi}{\omega_0})$.

Let

Clearly, for $\tau\in(\tau_k-\sigma_k, \tau_k+\sigma_k)$ and $(r, P)\in\partial\Omega_{\varepsilon_k, p_0}$ such that $\text{det} \Delta_{(E_+, \tau, p)}\Big(r+\text{i}\frac{2\pi}{p}\Big) = 0$ if and only if $\tau = \tau_k, p = p_0, r = 0$. This is the condition $(\text{A}_4)$.

So far, we have verified the conditions $(\text{A}_1)-(\text{A}_4)$ in ^[28]. Define

The condition ($\text{A}_4$) and $J(E_+, \tau_k, \frac{2\pi}{\omega_0}) = 1$ imply $H_1^\pm(E_+, \tau_k, \frac{2\pi}{\omega_0})(r, p)\neq0$ for $(r, p)\in\partial\Omega_{\varepsilon_k, p_0}$. Therefore, the first crossing number $\gamma_1(E_+, \tau_k, \frac{2\pi}{\omega_0})$ is calculated as

By the similar arguments, we can know $(E_0, \tau, p)$ and $(E_1, \tau, p)$ are not the centers, but $(E_2, \tau_k^*, \frac{2\pi}{\omega_0^*})$ and $(E_-, \tau_k, \frac{2\pi}{\omega_0})$ are isolated centers. we may also obtain $\gamma_1(E_2, \tau_k^*, \frac{2\pi}{\omega_0^*}) = -1$ and $\gamma_1(E_-, \tau_k, \frac{2\pi}{\omega_0}) = -1$.

In what follows we define

$\Sigma = Cl\{(z, \tau, p)|z \ is \ p-periodic \ solution \ of \ (5.1)\},$

$N = \{(\bar{z}, \tau, p)|F(\bar{z}, \tau, p) = 0\}.$

Let $C\Big(E_+, \tau_k, \frac{2\pi}{\omega_0}\Big)$ denote the connected component through $(E_+, \tau_k, \frac{2\pi}{\omega_0})$ in $\Sigma$.

By Theorem 3.2 in ^[28], there exists an isolated center $(E_+, \tau_k, \frac{2\pi}{\omega_0})$ which satisfies $J(E_+, \tau_k, \frac{2\pi}{\omega_0}) = 1$ and $\gamma_1(E_+, \tau_k, \frac{2\pi}{\omega_0})\neq0$ such that $C\Big(E_+, \tau_k, \frac{2\pi}{\omega_0}\Big)$ through $(E_+, \tau_k, \frac{2\pi}{\omega_0})$ in $\Sigma$ is nonempty. In addition, all the centers of (5.1) are isolated centers and satisfy

By Theorem 3.3 in ^[28], we obtain the following Lemma.

Lemma 5.1 $C\Big(E_+, \tau_k, \frac{2\pi}{\omega_0}\Big)$ is unbounded.

On the other hand, from Theorem 2.2 in Section 2, it is easy to obtain the following lemma.

Lemma 5.2 If $\tau$ is bounded, then any nontrivial periodic solution of system (1.2) is uniformly bounded.

Lemma 5.3 When the conditions (H) and $s < \delta$ hold, (1.2) has no any nontrivial $\tau-$periodic solutions.

Proof. If $(u^*(t), v^*(t))$ is a nontrivial $\tau-$periodic solution of (1.2), then it is also a nontrivial periodic solution of the following (5.5).

When the condition (H) hold, (5.5) has at most five equilibria, namely, $E_0 = (0, 0)$, $E_1 = (1, 0)$, $E_2 = (0, n\alpha\beta+\beta)$ and $E_\pm = (\bar{u}_\pm, \bar{v}_\pm)$. Note that $E_0$, $E_1$, and $E_2$ are located in $u$-axis and $v$-axis, they can not produce any nontrivial periodic solution. On the other hand, $E_\pm$ also can not produce any nontrivial periodic solution due to $s < \delta$. Thus, there is no any nontrivial periodic solution in (5.5). The proof is complete.

Theorem 5.4 Suppose the conditions (H), $G+D > 0$ and $s < \delta$ hold, then for each $\tau > \tau_k, k = 1, 2, 3\cdots$, system (1.2) has at least $k-1$ periodic solutions.

Proof. Let $\text{Proj}_\Theta C(E_+, \tau_k, \frac{2\pi}{\omega_0})$ be the projection of $C(E_+, \tau_k, \frac{2\pi}{\omega_0})$ onto $\Theta$-space. From Lemma 5.2, it is easy to know that $\text{Proj}_z C(E_+, \tau_k, \frac{2\pi}{\omega_0})$ is bounded. It follows from the proof of Lemma 5.3 that (1.2) with $\tau = 0$ has no nontrivial periodic solution. Thus, $\text{Proj}_\tau C(E_+, \tau_k, \frac{2\pi}{\omega_0})$ is away from zero.

We suppose that $\text{Proj}_\tau C(E_+, \tau_k, \frac{2\pi}{\omega_0})$ is bounded. Then there exist $\tau^*$ such that $\text{Proj}_\tau C(E_+, \tau_k, \frac{2\pi}{\omega_0})$ is a subset of interval $(0, \tau^*)$. From the definition of $\tau_k$ and $\omega_0$ in section 3, we have $\frac{2\pi}{\omega_0} < \tau_k, k = 1, 2, 3\cdots$. Applying Lemma 5.3, we have $p\in (0, \tau^*)$ when $(z, \tau, p)\in C(E_+, \tau_k, \frac{2\pi}{\omega_0})$.

The above discussion shows that $C(E_+, \tau_k, \frac{2\pi}{\omega_0})$ is bounded, which contradicts Lemma 5.1. Therefore, $\text{Proj}_\tau C(E_+, \tau_k, \frac{2\pi}{\omega_0})$ contains at least an interval $(\tau_k, +\infty)$. The proof is complete.

6.   Numerical simulations

Based on the previous discussion, three numerical results of system (1.2) are presented.

Case 1. We consider system (1.2) with $s = 0.8$, $m = 0.2$, $\alpha = 1$, $\beta = 0.6$, $\delta = 0.4$, $n = 0.8$, that is

We have $A = 0.4400$, $B = 0.0640$ and $\sqrt{A^2-4B} = -0.624$. It implies system (6.1) has no positive equilibrium and satisfies the condition (ii) in Theorem 3.2. Therefore, the prey-extinction equilibrium $E_2 = (0, 1.08)$ is stable when $\tau < \tau_0^*$, where $\tau_0^* = 2.1817$. When $\tau$ pass through the critical value $\tau_0^*$, $E_2$ loses its stability(see Figure 1).

Case 2. we discuss the following system

where $s = 0.2$, $m = 0.1$, $\alpha = 1$, $\beta = 0.6$, $\delta = 0.3$, $n = 0.4$. Calculate directly, we have $B = -0.5320$, $G = 0.7641$, $D = 0.2444$, then there is a unique coexisting equilibrium $E_+ = (0.8476, 1.1791)$. From the condition (ii) in Theorem 3.5, there exists the critical value $\tau_0 = 6.0476$ such that $E_+$ is stable when $\tau < \tau_0$ and unstable when $\tau > \tau_0$(see Figure 2).

On the other hand, we obtain the following values by using Matlab

$\alpha_1 = -6.9759 - 2.3462i$, $\alpha_2 = 2.7585 + 7.7731i$, $M = -0.0090 - 0.0008i$,

$g_{20} = -14.4293 -14.6239i$, $g_{02} = 19.8351 - 5.9477i$, $g_{11} = 0.6641 - 2.1798i$,

$g_{21} = -41.2805 +50.5284i$, $c_1(0) = -27.6348 -37.3972i$, $\lambda'(\tau_0) = 0.0192 - 0.0298i$.

It follows that $\mu_2 = 1440.4 > 0$ and $\beta_2 = -55.2696 < 0$ and $T_2 = 312.4076$, which, together with Theorem 4.1, implies that the bifurcating periodic solution exists when $\tau > \tau_0$ and the bifurcating periodic solution is stable on the center manifold and the period increases.

When $\tau = \tau_0$, all roots of characteristic equation (3.3) have negative real parts except $\pm i\omega_0$. Since the periodic solution on the center manifold is stabile, the periodic solution in the whole phase space is stable(see Figure 2 (c) and (d)).

Case 3. the simulation of the system (1.2) with $s = 0.6$, $m = 0.2$, $\alpha = 0.3$, $\beta = 0.5$, $\delta = 0.65$, $n = 0.6$ is given by

which has two positive equilibriums $E_- = (0.0145, 0.5987)$ and $E_+ = (0.2755, 0.7553)$ due to $A = -0.2900$, $B = 0.0040$ and $\sqrt{A^2-4B} = 0.2610$. By the condition (i) and (ii) in Theorem 3.5, we can know that $E_- = (0.0145, 0.5987)$ is always unstable for any $\tau > 0$ and there exists the critical value $\tau_0 = 1.2379$ such that $E_+$ is stable for any $\tau\in [0, \tau_0)$. When $\tau$ crosses $\tau_0$, $E_+$ is unstable and a Hop bifurcation occurs. We obtain $c_1(0) = -23.7642 +34.8292i$ and $\lambda'(\tau_0) = 0.5009 - 0.6450i$ by using Matlab, then $\mu_2 = 47.4457 > 0$ and $\beta_2 = -47.5284 < 0$ and $T_2 = -3.6496$, which implies that the bifurcating periodic solution is stable on the center manifold and the period increases. As discussed in system (6.2), for $\tau > \tau_0$, the periodic solution of system (6.3) in the whole phase space is stable. The corresponding simulation results are shown in Figure 3.

In addition, system (6.3) also has a prey-extinction equilibrium $E_2 = (0, 0.59)$ corresponding to the critical value $\tau_0^* = 0.6144$. Thus, $E_2 = (0, 0.59)$ and $E_+ = (0.2755, 0.7553)$ are the stability equilibria of system (6.3) when $\tau < 0.6144$, which is depicted in Figure 4.

7.   Conclusions

In the paper, we investigate a delayed predator-prey system with additional food and asymmetric functional response. The local stability of all possible equilibria are studied. It shows that we can exterminate the prey by adjusting the quality and quantity of additional food when the prey density and time delay are relatively small. We know that the number of positive equilibria is determined by the value of $B$. For $B < 0$, there is only one conditionally stable or unstable coexisting equilibrium $E_+$, which depends on the delay. However, there exist an absolutely unstable coexisting equilibrium$E_-$ and a conditionally stable or unstable coexisting equilibrium $E_+$ for $B > 0$ in system (1.2). We also find that the model is bistable when $B > 0, \ A < 0$ and $A^2-4B > 0$ (see Figure 4).

Our investigation shows that coexisting equilibrium $E_+$ is always unstable after $\tau$ passes through the first critical value $\tau_0$. That is to say, there does not exist any stability switching. However, stability switching can occurs in some predator-prey systems (see e.g. ^[26,29,33]). Moreover, the formulas determining the direction ($\mu_2$) and stability ($\beta_2$) of Hopf bifurcation are given. We also show that the local Hopf bifurcation implies the global Hopf bifurcation of positive equilibrium after the second critical value of delay. Finally, we give three examples to illustrate the stability of the system (1.2) near the first critical value, and the simulation results are consistent with our conclusions.

This paper mainly considers the effects of providing additional food and designing a delayed feedback on Holling-Tanner model theoretically. We cannot claim that our method always holds for Holling-Tanner model due to the variety of ways to provide additional food and design delayed feedback. However, our investigation has potential significance for biological control. Therefore, the future works may consider how does the different ways of additional food and delayed feedback affect a predator-prey system, and cover the effects of additional food on an ecoepidemic model and so forth.

Acknowledgments

The authors thank the editor and referees for their valuable suggestions and comments, which improved the presentation of this manuscript.

Conflict of interests

For the publication of this article, no conflict of interests among the authors is disclosed.