Global Hopf bifurcation of a delayed phytoplankton-zooplankton system considering toxin producing effect and delay dependent coefficient

1.   Introduction

Almost all aquatic life is based on plankton, the most abundant form of life, floating freely on well known aquatic surfaces such as wells, lakes, rivers, estuaries and oceans. The studies on plankton materials is an important area for the researchers who engaged in the field of marine ecology. The phytoplankton, the plant forms of plankton enables to engage with the process of photo-synthesis when sunlight placed on them and serves as the fundamental food source. The zooplankton, the animal forms of plankton consumes phytoplankton which is the most favorite food source of aquatic animals such as fish, shellfish, molluscs and jellyfish. One of the merits of the phytoplankton is that it can release oxygen by absorbing carbon dioxide, this process contributes to make a cleaner environment. On the other side, some phytoplankton can release toxin substances which are hazards for the aquatic animals, and even can kill the aquatic animals.

Under certain conditions, phytoplankton can develop rapidly. These bodies of water which are enriched in phytoplankton may become red tide, while the predator of zooplankton will have a certain impact on the number of phytoplankton, so it is of great significance to study their relationship for the occurrence of phytoplankton bloom. Mathematical models have made considerable contributions to the understandings of the relationship between phytoplankton and zooplankton [1,2,3,4,5,6,7,8,9,10,11,12,13]. Toxin available in water sources impact on plankton materials and also can be used to suppress the growth of plankton blooms. J. Chattopadhyay et al. [14] proposed a mathematical model to study the behavior of toxic-phytoplankton (TPP) and zooplankton as well as study their interactions. The general form of mathematical model proposed in [14] can be governed by the following two nonlinear ordinary differential equations,

where $X(t)$ represents the density of TPP population and $Y(t)$ is the density of zooplankton population at time $t$. In system (1.1), $r$ represents the growth rate, and $A$ is the environmental carrying capacity for TPP. The functional response functions for zooplankton grazing phytoplankton and the toxin distribution causing zooplankton death are denoted by $\psi(X)$ and $\phi(X),$ respectively. Further, $\beta$ is the maximum uptake rate and $\beta_1$ is the ratio of conversion satisfying $\beta_1 < \beta$. The natural death rate of zooplankton is $\nu$. The parameter $\rho$ represents the rate of toxin release of phytoplankton. All parameters are positive.

We know that the delay caused by the maturity of TPP plays a crucial role on the dynamic behavior of phytoplankton zooplankton system, which seems that delay could cause rich dynamics [15,16,18,17,19]. In phytoplankton zooplankton system, many researchers assumed that the releasing of toxin is an instantaneous process, while the real case is that time delay may have an essential influence on the dynamics of mathematical models. In [16], Chattopadhyay et al. analyzed the following delayed system

In system (1.2), functions $\psi$ and $\phi$ are assumed to $\psi(X) = X$ and $\phi(X) = X/(\gamma + X).$ The process of toxin releasing follows a discrete variation and $\tau$ is the maturation time of phytoplankton for releasing toxin, the authors studied the oscillation behavior of phytoplankton and zooplankton populations based on system (1.2) with a condition that helps the oscillatory behavior to be stable. By regarding delay as the bifurcation parameter, Saha and Bandyopadhyay [20] modified system (1.2) by using Holling type Ⅱ function instead of $\psi(X)$ and by discussed the oscillatory behavior of phytoplankton and zooplankton. Moreover, the authors also established a set of conditions for the existence of globally periodic solutions. In the view of real ecological meaning, a gestation period of zooplankton is likely to delay the contact with phytoplankton. Rehim and Imran [21] put on a more general model:

where the first term and the final term with delay in the second equation of system define the zooplankton's gestation delay and TPP's maturity delay, respectively. The authors obtained the globally asymptotical stability properties of equilibria and established some conditions such that occurring of local Hopf bifurcation at the positive equilibrium. Based on Rehim and Imran [21], Wang et al. [22] introduced the harvesting term in zooplankton population, and the globally asymptotical stability properties of equilibrium and occurrence of Hopf bifurcation are obtained. They also established excellent result that the system had at least one positive periodic solution when delay changes. Furthermore, Jiang et al. [23] replaced the functions $\psi(X)$ and $\phi(X)$ by Holling Ⅱ type functions with the different half-saturation constants. They established the globally asymptotical stability properties of boundary equilibrium and the existence of local and global Hopf bifurcation under certain conditions.

Motivated by the systems considered in [21,22,23] and in this paper, we formulated a generalized system as shown below,

where $\phi(X)$ conforms to the following hypothesis:

The prototypes of response function $\phi(X)$ can be found in the literatures[24,25,26,27], for example, $\phi(X) = X$ (Holling type Ⅰ or Lotka-Volterra kinetics), $\phi(X) = X/(m+X)$ (Holling type Ⅱ or Michaelis-Menten kinetics), $\phi(X) = X^2/(m+X^2)$ (Holling type Ⅲ) and, $\phi(X) = 1-e^{-\alpha X}$ (Ivlev's functional response).

This paper is organized as follows. In section 2, it will give the positivity and boundedness of solutions of system (1.4). In section 3, a special case of system (1.4) that is without considering time delay will be properly studied and also the stability properties of each of the equilibria have been extensively presented. In section 4, we obtain results related to stability of equilibria and existence of Hopf bifurcation of system (1.4) with delay. In section 5, employing the center manifold theory that are presented in scholar manuscript [28,29,30], the properties of Hopf bifurcation are derived. For delayed systems, are there exist of large-scale periodic solutions when $\tau$ increases from the first Hopf bifurcation values? Authors [31,32] studied the global Hopf bifurcation of a delayed system by using a the result due to Wu [33]. In the following, by using the global Hopf bifurcation result [33], we obtain the results that are related the global existence of periodic solution with delay varying in section 6. In section 7, numerical simulations are given with the purpose of verifying the theoretical results. At last, some conclusions are given in section 8.

2.   Positivity and boundedness of solutions

Let $n$ be an integer, for any $n\geq 1$, define $\mathbf{R}_+^n = \{(x_1, x_2, \cdots, x_n)\in \mathbf{R}^n, x_i\geq 0, 1\leq i \leq n\}$ and $\mathbf{Int}\; \mathbf{R}_+^n = \{(x_1, x_2, \cdots, x_n)\in \mathbf{R}^n, x_i > 0, 1\leq i \leq n\}$. For $\tau > 0$, denote the space $C([-\tau, 0], \mathbf{R}^2)$ with the supermum norm. Let $C_+^2 = C([-\tau, 0], \mathbf{R}_+^2)$ and $\mathbf{Int}\; C_+^2 = C([-\tau, 0], \mathbf{Int\; R}_+^2).$ For system (1.4), let us take any initial function on $\mathbf{Int}\ C_+^2$ and it always assume that $\beta > \beta_1.$

Lemma 2.1. Let $X(\theta): = \varphi_1(\theta), Y(\theta): = \varphi_2(\theta)\geq 0$ for $\theta\in[-\tau, 0)$, and $\varphi_1(0), \varphi_2(0) > 0$. Then for some $\sigma > 0$, all solutions of system (1.4) uniquely exist on $[0, \sigma)$ when the initial function $\varphi = (\varphi_1, \varphi_2)^T\in \mathbf{Int}\ C_+^2$, and all solutions are positive for $t\in [0, \sigma)$.

Proof. By [34], solutions of system (1.4) with initial function $\varphi\in \mathbf{Int}\; C_+^2$ uniquely exist on $t\in [0, \sigma)$ for some $\sigma > 0$. Suppose that $(X(t), Y(t))$ is the solution of (1.4) for $t\in [0, \sigma)$. Without loss of generality, it assumes that $[0, \sigma)$ is the maximum existence internal of solution. If the solution exists for all $t > 0$, then $\sigma = \infty$. Integrating the first equation of system (1.4) gives

Next, the proof by contradiction is used to show $Y(t) > 0$. Suppose that there exists $t^*\in[0, \sigma)$ such that

In the equation of $Y(t)$, exchanging $t$ for $t^*$, it has

It is a contradiction to $\dot{Y}(t^*)\leq 0$. So that, $Y(t) > 0$ for all $t\in [0, \sigma).$ This completes the proof.

Lemma 2.2. The nonnegative solutions of system (1.4) are bounded on the interval $t\in [0, \sigma)$. Further, $\lim\limits_{t\rightarrow\infty} \sup X(t)\leq A$.

Proof. Let $(X(\theta), Y(\theta))$ be the solution of system (1.4). From the equation of $X(t)$ in system (1.4), it has that $\dot{X}(t)\leq \dfrac{r}{A}X(t)(A-X(t))$, which yields that $\lim\limits_{t\rightarrow\infty} \sup X(t)\leq A$. Define

Then

By the comparison theory [35], $\mathcal{W}(t)\leq \digamma(t)$, where $\digamma(t) = \digamma(0)e^{-\nu t}+\dfrac{\beta_1A}{4\beta r\nu }e^{-\nu \tau}(\nu +r)^2(1-e^{-\nu t})$ is the solution of initial value problem

Consequently, $\mathcal{W}(t)\leq \mathcal{W}(0)+\dfrac{\beta_1 A}{4\beta r\nu }e^{-\nu t}(\nu +r)^2$. Furthermore, it has

Up to now, it has obtained that the nonnegative solutions of system (1.4) are bounded on the interval $t\in [0, \sigma).$ This completes the proof.

According to the continuation theorem of solutions for FDE [34], we can get the following theorem.

Theorem 2.1. The solutions $(X(t), Y(t))$ of system (1.4) with the initial function $\varphi\in \mathbf{Int}\; C_+^2$ is uniquely existent, nonnegative and bounded on $[0, +\infty)$ and satisfies

where

The following result is given for our need.

Lemma 2.3. For the following system

if $\mathfrak{a} < \mathfrak{b}$, then the solution satisfies $\lim\limits_{t\rightarrow\infty}\mathfrak{y}(t) = 0,$ where $\mathfrak{a, b}, \tau > 0$ and $\mathfrak{y}(t) > 0$ for $t\in [-\tau, 0]$.

3.   The stability analysis of system(1.4) without delay

In this part, it investigates the stability of ODE system, that is, $\tau = 0$. When $\tau = 0$, system (1.4) has three equilibria $E_0(0, 0), E_1(A, 0)$ and $E^*(X^*(0), Y^*(0))$, where $X^*(0)$ and $Y^*(0)$ satisfy $r(A-X^*(0))-A\beta Y^*(0) = 0,$ and $\beta_1 X^*(0)-\rho \phi (X^*(0)) = \nu$. The positive equilibrium $E^*$ is existent and unique if $\beta_1 A-\rho \phi (A) > \nu$. We can obtain $E_0$ is always a saddle point. If $\beta_1 A-\rho \phi (A) < \nu$, then $E_1$ is locally asymptotically stable (LAS), and, if the inequality is reversed, then $E_1(A, 0)$ is unstable and $E^*$ is locally asymptotically stable. Furthermore, it has the global stability result of $E_1$.

Theorem 3.1. If $\beta_1 A-\rho \phi (A) < \nu$, then $E_1$ is globally asymptotically stable (GAS).

Proof. By the positivity of solutions, it has that $\dot{X}(t)\leq rX(t)(1-X(t)/A)$. This implies that $\lim\limits_{t\rightarrow\infty}\sup X(t) \leq A.$ Consequently, $\dot{Y}(t) = -\nu Y(t)+\beta_1X(t)Y(t)-\rho \phi (X(t))Y(t) = [\nu +\beta_1 X(t)-\rho \phi (X)]Y(t)$ $\leq$ $[-\nu +\beta_1 A -\rho \phi (A)]Y(t)$. Let $y(t)$ be the solution of $\dot{\mathbf{y}}(t) = [-\nu +\beta_1 A-\rho \phi(A)]\mathbf{y}(t)$. Since $\beta_1 A-\rho \phi (A) < \nu$, it has $\mathbf{y}(t)\rightarrow 0$. Hence, $Y(t)\rightarrow 0$. By comparison theorem, $Y(t)$ is bounded. Let $\eta\in(0, A)$, then there exists $T_\eta > 0$ such that $\dot{X}(t)\geq \dfrac{r}{A}X(t)(A-\eta-X(t))$ for $t\geq T_\eta$, hence, we have $\lim\limits_{t\rightarrow\infty}\inf X(t)\geq A-\eta.$ Since $\eta\in(0, A)$ is arbitrary, $\lim\limits_{t\rightarrow\infty}\inf X(t)\geq A.$ We already have $\lim\limits_{t\rightarrow\infty}X(t)\leq A.$ Then, $\lim\limits_{t\rightarrow\infty}X(t) = A.$ Hence $E_1$ is GAS.

4.   The dynamic analysis of system (1.4) with delay

Now let us return the case $\tau\geq 0$. System (1.4) has three equilibria $E_0(0, 0), E_1(A, 0)$ and $E^*(X^*(\tau), Y^*(\tau))$, where $X^*(\tau), Y^*(\tau)$ satisfy $r(1-X^*(\tau)/A)-\beta Y^*(\tau) = 0, \beta_1 X^*(\tau)-\rho \phi (X^*(\tau)) = \nu e^{\nu \tau}$. Therefore, we have that $E^*$ exists and is unique if $0\leq\tau < \tau_c\doteq(1/\nu)\ln[(1/\nu)(\beta_1 A-\rho \phi (A))]$ and $E^*\neq E_1.$ If $\tau > \tau_c$, then $E^*$ is not existent.

That is clear that $E_0$ is always a saddle point. To $E_1,$ the characteristic equation of system (1.4) at $E_1$ is given by

Then we have the following theorem.

Theorem 4.1. If $0\leq\tau < \tau_c$, then $E_1$ is unstable. If $\tau > \tau_c,$ then $E_1$ is GAS.

Proof. One of the roots of (4.1) is $\lambda = -r,$ and the other roots satisfy

It has $\mathbf{H}(0) = \nu e^{\nu \tau}, \dot{\mathbf{H}}(\lambda) > 0, \mathbf{H}(+\infty) = \infty$. Since $\tau < \tau_c$, there exists a unique positive root $\lambda$ such that (4.2) holds and (4.1) has at least one positive root $\lambda$. Hence, $E_1$ is unstable when $\tau\in [0, \tau_c).$

Now, one proves $E_1$ is GAS when $\tau > \tau_c$. From the second equation, it has that $\dot{Y}(t) = -\nu Y(t)+e^{-\nu \tau}[\beta_1X(t-\tau)-\rho \phi (X(t-\tau))]Y(t-\tau)\leq -\nu Y(t)+e^{-\nu \tau}[\beta_1 A-\rho \phi (A)]Y(t-\tau)$. Let $\mathbf{y}(t)$ be the solution of $\dot{\mathbf{y}}(t) = -\nu \mathbf{y}(t)+e^{-\nu \tau}[\beta_1 A-\rho \phi (A)]\mathbf{y}(t-\tau)$. Since $\tau > \tau_c,$ Lemma 3 guarantees that $\mathbf{y}(t)\rightarrow 0$. Hence, $Y(t)\rightarrow 0$. Furthermore, it can get that $X(t)\rightarrow A,$ which implies the global stability of $E_1.$ This completes the proof.

Let

then

Let $\lambda = \pm i\omega^0(\tau)\text{ } (\omega^0(\tau) > 0)$ be a pair of simple pure imaginary roots of (4.1). The stability switches may occur at the $\tau$ values:

where

and $\theta^0_{\pm}(\tau)\in[0, 2\pi)$ satisfy

where $R = e^{-\nu \tau}[\beta_1A-\rho \phi (A)]$. Hence, $\omega_-^0(\tau)$ is infeasible; $\omega_+^0(\tau)$ is feasible if $\tau\in [0, \tau_c)$ and $\tau_n^+(\tau) = (1/\sqrt{R^2-\nu ^2})\Big\{2\pi-\text{arcos}(\nu /R)+2n\pi\Big\}$; $\omega_+^0(\tau)$ is infeasible for $\tau > \max\{0, \tau_c\}$.

Define

Theorem 4.2. If $0\leq\tau < \tau_c$, then $\pm i\omega_+^0(\tau)$ are a pair of simple pure imaginary roots of system (1.4), and $Z_n(\tau_*) = 0$ for some $n\in \mathbb{N}$ and some $\tau_*\in [0, \tau_c)$. This pair of roots cross the imaginary axis from left (right) to right (left) if $\kappa_+^0(\tau_*) > 0\; (<0)$, where

If $\mathbb{J}^0$ is not empty. For all $\tau_j\in \mathbb{J}^0$, if $Z_{n}'(\tau_j)\neq 0$ holds, then system (1.4) undergoes Hopf bifurcations at $E_1$ when $\tau = \tau_j$.

Remark 4.1. It has proved that there are only two equilibria $E_0$ (unstable) and $E_1$ (GAS) when $\tau > \tau_c$. However if $\tau\in [0, \tau_c)$ and $\beta_1A-\rho \phi (A) > \nu$, then both $E_0$ and $E_1$ are unstable, at the same time, $E^*$ exists and is LAS.

For convenience, we denote $X^*(\tau) = X^*, Y^*(\tau) = Y^*$. Let $x = X-X^*, y = Y-Y^*$, then system (1.4) becomes

whose characteristic equation is

where

Equation (4.4) converts to the general form

where

Equation (4.5) takes all the coefficients of $\mathbb{P}$ and $\mathbb{Q}$ depending on $\tau.$ We apply the criterion due to Bereta and Kuang [36], let $\lambda = i \omega\; (\omega = \omega(\tau) > 0)$ be a root of (4.5), then $\omega$ satisfies

and

By the definitions of $\mathbb{P}$ and $\mathbb{Q},$ and applying the property (i), (4.6) becomes:

which yields

and its roots are given by

where

Let

If $I$ is non-empty, then $\omega = \omega_+$ and $\omega_-$ is not feasible for all $\tau\in I$, and $\omega$ is not definited for $\tau\notin I.$

Then, for $\tau\in I$, let $\theta(\tau)\in [0, 2\pi)$ be defined by

Hence, it can define the maps $\tau_{n}(\tau)$ given by

Next, let us introduce the following continuous and differentiable functions

Theorem 4.3. Let $\lambda = \pm i \omega_+(\tau^*)$ be a pair of simple pure imaginary roots of Eq (4.5), and at some $\tau^*\in I$,

This pair of roots crosses the imaginary axis from left (right) to right (left) if $k_+(\tau^*) > 0\; (<0)$, where

Remark 4.2. For $\tau\in I$, it has the following one sequence of functions:

Clearly $S_{n}^+(\tau) > S_{n+1}^+(\tau)$ for all $n\in \mathbb{N}, \; \tau\in I.$

Theorem 4.4. For system (1.4), if $\beta_1A-\rho \phi (A) > \nu$, then it has the following conclusions.

(i) If $I = \emptyset$ or non-empty set, but $S_{n}(\tau) = 0$ has no positive root in $I$, then $E^*$ is LAS for all $\tau\in [0, \tau_c)$;

(ii) If $I\neq \emptyset$ and $S_{n} = 0$ has positive roots in $I$, denoted by $\tau_n^j$, for some $n\in \mathbb{N}$. it assumes that $S_n'(\tau_n^j)\neq 0$. Rearrange these roots as the set $\mathbb{J}: = \{\tau^0, \tau^1, ..., \tau^m\}$ with $\tau^j < \tau^{j+1}, j = 0, ..., m-1.$ Then $E^*$ is LAS for $\tau\in[0, \tau^0)\cup (\tau^m, \bar{\tau})$ and unstable for $\tau\in (\tau^0, \tau^m).$ Hopf bifurcations occur at $E^*$ when $\tau = \tau^j$.

5.   The direction and stability of Hopf bifurcation at $E^*$

In section 4, we have obtained the sufficient conditions which ensure that system (1.4) undergoes Hopf bifurcation at $E^*.$ In this section, using the center manifold theory presented by Hassard et al. [28], we can establish an explicit formula to determine the bifurcating direction and stability of periodic solutions at $\tau = \tau^0$. In fact, let $\omega^0 = \omega(\tau^0)$. Furthermore, let $\tau = \tau^0+\mu$, then $\mu = 0$ is a Hopf bifurcation value of system (1.4). We rewrite (4.3) as

where

and for $\varphi\in C$, define

By the Riesz representation theorem, there exists a matrix $\zeta(\theta, \mu)$ in $\theta\in[-\tau, 0]$ whose elements are bounded variation functions such that

In fact, it can choose

where

Choose $\varphi\in C^1 = C([0, \tau], (\mathbf{R}^2)^{*})$, define

and

where

Let $\mathfrak{u} = (x, y)^T,$ then system (5.1) becomes

For $\psi\in C^1$, define

and a bilinear form

Then $\mathbb{A}^*$ and $\mathbb{A}$ are adjoint operators, and they have the same eigenvalues $\pm i\omega^0$. By computation, we obtain that the eigenvector of $\mathbb{A}(0)$ corresponding to $i\omega^0$ is $\mathfrak{q}(\theta) = \Big(1, \mathfrak{q}_2\Big)e^{i\omega^0\theta}$ and the eigenvector of $\mathbb{A}^*$ corresponding to $-i\omega^0$ is $\mathfrak{q}^*(s) = \bar{D}\Big(1, \mathfrak{q}_2^*\Big)^Te^{i\omega^0 s}$, where

Moreover,

where

Let $\mathfrak{ u}_t$ be the solution of (5.1) with $\tau = \tau^0$. Define $\mathfrak{z}(t) = \langle \mathfrak{q}^*, \mathfrak{u}_t \rangle, \mathfrak{u}_t = (x_t, y_t)$, then

where

Rewriting (5.3) as

where

Substituting (5.2) and (5.3) into $\dot{W} = \dot{\mathfrak{u}}_t-\dot{\mathfrak{z}}\mathfrak{q}-\dot{\bar{\mathfrak{z}}}\mathfrak{q}$, it has

where

By comparing the coefficients, it yields

For

then

Notice that

Comparing the coefficients, it can obtain

It still needs to compute $W_{20}(\theta)$ and $W_{11}(\theta).$ When $\theta\in [-\tau, 0),$

By comparing the coefficients with (5.4), it can yield that

In addition,

Furthermore,

and similarly

where $\mathcal{E}_1$ can be determined by setting $\theta = 0$ in $\mathbb{H}$. In fact,

where

Hence,

and

Notice that

hence,

Similarly, it has

Furthermore, we get

and

Then $g_{21}$ can be computed. Hence each $g_{ij}$ can be determined by the parameters. Thus, the following quantities can be computed:

Hence, we have the following theorem.

Theorem 5.1. If $\mathcal{U}_2 > 0\; (<0)$, the direction of Hopf bifurcation is supercritical (subcritical); if $\mathcal{B}_2 < 0\; (>0)$, the bifurcation periodic solutions are orbitally stable (unstable); if $\mathcal{T}_2 > 0\; (<0)$, the period increase (decrease).

6.   Global existence of Hopf bifurcations

In this part, it will consider the global existence of periodic solution. Firstly, one gives some notations [37].

Define

It assumes that $\mathbb{J}_+\neq \emptyset$. In the following, the global existence of periodic solutions bifurcating from $(E^*, \tau^j, \frac{2\pi}{\omega_j^+})$ for system (1.4) is investigated by using global Hopf bifurcation theorem [33], where $\omega_j^+ = \omega^+(\tau^j)$ ($\tau^j\in \mathbb{J}^+$), and $\pm i\omega_j^+$ is a pair of simple roots of (4.5) when $\tau = \tau^j$.

Next, it always assumes that $\tau\in [0, \tau_c)$ and $\beta_1 A-\rho \phi (A) > \nu$ are satisfied. For convenience, let $\mathfrak{Z}_t = (X_{t}, Y_{t})$, system (1.4) can be rewritten as the following FDE:

where $\mathfrak{Z}_t(\theta) = \mathfrak{Z}(t+\theta).$ By the results in [33], it can define the following signs:

and let $\mathcal{C}_{(\mathfrak{Z}^*, \tau^j, \frac{2\pi}{\omega_j^+})}$ denote the connected component of $(\mathfrak{Z}^*, \tau^j, \frac{2\pi}{\omega_j^+})$ in $\mathfrak{L}$, and $\text{Proj}_{\tau}(\mathfrak{Z}^*, \tau^j, \frac{2\pi}{\omega_j^+})$ represents its projection on $\tau$ space, where $\mathfrak{Z}^* = E^*.$ On the basis of Theorem 2.1, it has

Lemma 6.1. All periodic solutions of system (6.1) are uniformly bounded in $\mathbf{R}_+^2.$

Lemma 6.2. System (6.1) does not exist non-trivial $\tau$-periodic solution.

Proof. It can know that $X$-axis and $Y$-axis are the invariable manifold of system (6.1) and the trajectories of system (6.1) are disjoint one another, which means that if there are the periodic solutions in the first quadrant, then $E^*$ must be in its interior. Using the proof by contradiction,

It assumes that system (6.1) has non-trivial $\tau$-periodic solution in the first quadrant, then the following system has the non-trivial periodic solution:

For another, define Dulac function $\mathcal{Q} = 1/(XY)$, then

Therefore, according to Dulac theorem, system (6.2) has no any periodic solution in the first quadrant, which leads to a contradiction.

Theorem 6.1. If $\mathbb{J}_+\neq \emptyset$, $\tau\in [0, \tau_c)$ and $\beta_1 A-\rho \phi (A) > \nu$ hold, then for each $\tau^j\in \mathbb{J}_+$, there must be $\tau^i\in \mathbb{J}^0\cup \mathbb{J}-\{\tau^j\}$, such that system (6.1) exists at least one positive periodic solution when $\tau$ varies between $\tau^i$ and $\tau^j$.

Proof. Firstly, it will prove that $\text{Proj}_{\tau}(\mathfrak{Z}^*, \tau^j, \frac{2\pi}{\omega_j^+})$ is unbounded for each $j.$ The characteristic matrix of system (6.1) at $\bar{\mathfrak{Z}} = (\bar{\mathfrak{Z}}^{(1)}, \bar{\mathfrak{Z}}^{(2)})\in \mathbf{R}_+^2$ satisfies:

System (6.1) has three equilibria $\bar{\mathfrak{Z}}_1, \bar{\mathfrak{Z}}_2$ and $\mathfrak{Z}^*$. It knows from section 4 that $(\bar{\mathfrak{Z}}_1, \tau, \mathfrak{p})$ is not center while $(\bar{\mathfrak{Z}}_2, \tau, \mathfrak{p})$ and $(\mathfrak{Z}^*, \tau, \mathfrak{p})$ are isolated centers. Base on Theorem 4.4 and the condition $S_n'(\tau^j)\neq 0$, there exist $\epsilon > 0$, $\upsilon > 0$ and a smooth curve $\Lambda:$ $(\tau^j-\upsilon, \tau^j+\upsilon)\rightarrow C$, such that, $\det(\triangle (\Lambda(\tau^j))) = 0$, $|\Lambda(\tau^j)-i\omega_j^+ | < \epsilon$ for all $\tau\in[\tau^j -\upsilon, \tau^j+\upsilon]$ and $\Lambda(\tau^j) = i\omega_j^+$, $\mathit{\boldsymbol{dRe}}\Lambda(\tau^j)/{\rm{d}}\tau\neq 0.$

Let $\Omega_{\epsilon, \frac{2\pi}{\omega_j^+}} = \{(\eta, \mathfrak{p}):0 < \eta < \varepsilon, |\mathfrak{p}-\frac{2\pi}{\omega_j^+}| < \varepsilon\},$ then, on $[\tau^j-\upsilon, \tau^j+\upsilon]\times\partial\Omega_{\varepsilon, \frac{2\pi}{\omega_j^+}}$, $\Delta(\mathfrak{Z}^*, \tau, \mathfrak{p})(\eta+\frac{2\pi}{\mathfrak{p}}i) = 0$ iff $\eta = 0, \tau = \tau^j, \mathfrak{p} = \frac{2\pi}{\omega_j^+}$. Furthermore, define

then it has the crossing number:

Hence it has

where $(\bar{\mathfrak{Z}}, \tau, \mathfrak{p})$ takes the form of either $(\mathfrak{Z}^*, \tau^j, \frac{2\pi}{\omega_j^+})$ or $(\bar{\mathfrak{Z}}_2, \tau^j, \frac{2\pi}{\omega_j^+}).$ Hence the connected component $\mathcal{C}_{(\mathfrak{Z}^{*}, \tau^{j}, \frac{2\pi}{\omega_j^+})}$ through $(\mathfrak{Z}^*, \tau^j, \frac{2\pi}{\omega_j^+})$ in $\mathfrak{L}$ is unbounded [33].

By representation of $\tau^j$, there exists a $j > 0$ such that

Therefore, it has that $\frac{\tau}{j+1} < \mathbb{T} < \tau$ if $(\bar{\mathfrak{Z}}, \tau, \mathbb{T})\in \mathcal{C}_{(\mathfrak{Z}^*, \tau^j, \frac{2\pi}{\omega_j^+})}$, which shows that the projection of $\mathcal{C}_{(\mathfrak{Z}^*, \tau^j, \frac{2\pi}{\omega_j^+})}$ is bounded onto the $\mathbb{T}$-space. On the other hand, Lemmas 4 and 5 imply that the projection of $\mathcal{C}_{(\mathfrak{Z}^*, \tau^j, \frac{2\pi}{\omega_j^+})}$ is bounded for $\tau\in \text{Proj}_{\tau} (\mathfrak{Z}^*, \tau^j, \frac{2\pi}{\omega_j^+})\cap(A_j, B_j)$ onto the $\mathfrak{Z}$-space. Therefore, either $[A^j, \tau^j]\subset \text{Proj}_{\tau} (\mathfrak{Z}^*, \tau^j, \frac{2\pi}{\omega_j^+})\cap (A_j, B^j)$ or $[\tau^j, B^j]\subset \text{Proj}_{\tau} (\mathfrak{Z}^*, \tau^j, \frac{2\pi}{\omega_j^+})\cap(A_j, B_j)$ holds. If not, $\mathcal{C}_{(\mathfrak{Z}^*, \tau^j, \frac{2\pi}{\omega_j^+})}$ is unbounded and $\text{Proj}_{\tau} (\mathfrak{Z}^*, \tau^j, \frac{2\pi}{\omega_j^+})\subset(A_j, B_j),$ which leads to a contradict.

7.   Numerical simulations

In this section, we shall use Matlab to perform some numerical simulations on system (1.4). We choose $\phi (X) = X$ and the parameter values as follows:

By computing, we obtain $\tau_c\doteq 5.9915$. From Figures 1 and 2, we can see that, for $n = 0$, $S_n(\tau) = 0$ has two roots and $\theta(\tau)$ intersects $\tau\omega_+(\tau)$ twice at $\tau^0\doteq0.123$ and $\tau^1\doteq 3.5462$. When $n\geq 1$ and $\tau\in [0, \tau_c)$, $S_n(\tau) = 0$ has no roots and $\theta(\tau)+2n\pi$ has no point of intersection with $\tau\omega_+(\tau).$ We choose the initial functions $X(t) = 2$ and $Y(t) = 1$ for $t\in[-\tau, 0]$.

Under the parameters (7.1), one can easily verify that $E^*(2.5, 1.9)$ is stable in absence of delay. For a small delay ($\tau = 0.05$), the equilibrium $E^*$ is stable (see Figure 3). However, when $\tau$ increases to $\tau^0\doteq 0.123$, $E^*$ becomes unstable because a Hopf bifurcation occurs at the moment. At the same time, a periodic orbit surrounding $E^*$ produces (see Figure 4 for $\tau = 0.13$), and by computing in sections 4 and 5, we can obtain $\omega^0 = 0.9415, \mathbf{Re} C_1(0) = -273.9906 < 0$, hence, the bifurcation periodic solution is stable. When $\tau$ changes from 0.13, the periodic solution changes its shape slightly (see Figure 5). As $\tau$ increasing, it returns to a periodic solution (see Figure 6). When $\tau\doteq4$, the periodic orbit disappears and $E^*$ restores stability (see Figure 7) and remains stability until $\tau\doteq 5.9915$. When $\tau = 7 > 5.9915$, zooplankton population dies out (see Figure 8). It can find that these results are consistent with the ahead theoretical analyses.

8.   Conclusions

In this paper, a system of describing the relationship between TPP and zooplankton is investigated. The process of toxin liberation uses the general function $\phi (X)$ with the restrictions (1.5). We induce two discrete delays to the consume response function and distribution of toxin term to explore the dynamical behavior as delay varying. By analysing the ODE system, we obtains the parameters conditions for the asymptotical stability of equilibrium. At the same time, if the natural death rate of zooplankton exceeds the threshold parameter, then zooplankton will die out ultimately and phytoplankton will persist, which means that phytoplankton bloom may break out. On the contrary, if the natural death rate of zooplankton is less than the threshold parameter, then phytoplankton and zooplankton will persistent coexistence and the number remains at a certain level, which means phytoplankton bloom can't occur. In this situation, the natural death rate of zooplankton is an important factor to the occurrence of phytoplankton bloom.

By analysing the system with time delay, we find that system (1.4) may occur the stable switches as delay changing. And Theorem 4.4 shows that although system (1.4) undergoes stability switches, system (1.4) is ultimately stable, which is significantly different from the stability switching phenomenon of a time-delay system with coefficients without parameters. Moreover, the global existence of periodic solutions is obtained, that is, under certain conditions, multiple positive periodic solutions will exist as delay changing. The results show that the delay (gestation delay from zooplankton and maturity delay of TPP) is a key factor to the periodic outbreaks of phytoplankton bloom. If the delay is more than some value, then zooplankton will die out and phytoplankton bloom can occur. If the delay is less than some value, then under the certain conditions, phytoplankton number will have periodic oscillations, which means the periodic outbreaks of phytoplankton bloom. System (1.4) generalize the system in [21,22,23] and these results obtained in this paper can also apply to above systems.

Acknowledgments

Z. Jiang was supported by National Natural Science Foundation of China (Nos. 11801014 and 11875001), Natural Science Foundation of Hebei Province (No. A2018409004), Support Plan for Hundred Outstanding Innovative Talents of Colleges and Universities in Hebei Province and Talent Training Project of Hebei Province; T. Zhang was supported by SDUST Research Funds (No. 2014TDJH102) and Scientific Research Foundation of Shandong University of Science and Technology for Recruited Talents.

Conflict of interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.