Optimal harvesting of a competitive n-species stochastic model with delayed diffusions

1.   Introduction

Plankton includes plants and animals that float freely in some fresh water bodies, and almost all aquatic life is based on plankton ^[1]. Aquatic ecosystems are affected by many factors, including physical and chemical signals in the environment, plankton predation and competition ^[2,3]. Many scholars have carried out analysis concerning the impact of environment on the ecosystem and the treatment of sewage ^[4,5,6]. We also know how important plankton itself is to the wealth of marine ecosystems and ultimately to the planet itself. On the one hand, plankton species have positive effects on the environment, such as providing food for marine life, oxygen for animal life; on the other hand, it have harmful effects, such as economic losses to fisheries and tourism due to algae blooms ^[7,8].

In recent years, different models of plankton have been established and studied, for example, model with two harmful phytoplankton ^[9], models with time delays ^[10,11] and stochastic models ^[12,13,14]. Toxins produced by harmful phytoplankton tend to be concentrated at higher levels in the food web, as they can spread through the marine food web, affecting herbivores at higher nutrient levels, reaching fish, and through them eventually reaching marine mammals, even in seabirds ^[9,11]. There is also some evidence that the occurrence of toxin-producing phytoplankton is not necessarily harmful, but rather helps maintain a stable balance of nutrient dynamics through the coexistence of all species. These results suggest that toxin-producing phytoplankton (TPP) play an important role in the growth of zooplankton populations ^[15].

It was shown that aquatic plant systems not only have extraordinary memories of climatic events, but also exhibit phenomenological responses based on memory ^[16,17,18]. The authors noted that environmental factors often alter the expression of chromatin in multiple responsive genes in ^[19]. Environment-induced chromatin markers are at certain sites and is transmitted by cell division, allowing plants to acquire memories of environmental experiences. This ensures that the plant can adapt to changes in its environment or perform better if the event occurs again. In some cases, it is passed on to the next generation, namely, epigenetic mechanisms. This mechanism is crucial for plants' stress memory and adaptation to the environment, suggesting that plants do form memory and defense mechanisms in certain environments. In addition, a large amount of zooplankton chemical signal learning and corresponding reactions have been documented for aquatic systems ^[20,21]. In summary, such memory and genetic characteristics can not be neglected for plankton systems.

As we all know, fractional order derivatives are a good tool for describing the memory and genetic properties of various materials and processes. In other words, the application of fractional order dynamical systems can fully reflect some long-term memory and non-local effects. That is, fractional differential equations have an advantage over classical integer differential equations for describing such systems. In recent years, more and more researchers began to study the qualitative theory and numerical solution of fractional order biological model ^[22,23,24]. The main reason is that fractional order equations are naturally related to memory systems that exist in most biological systems ^[25,26]. In addition, fractional-order derivative has also been widely studied and applied in physics^[27], engineering ^[28], biology ^[29] and many other fields ^[30,31,32]. At present, there are more than six definitions of fractional derivative, among which Riemann-Liouville and Caputo derivatives are the most commonly used ^[33]. In the case of fractional Caputo derivative, the initial conditions are expressed by the values of the unknown function and its integer derivative with clear physical meaning ^[34]. So we will adapt the Caputo's definition in our paper.

The interactions between phytoplankton and zooplankton do not occur instantaneously in real ecosystems. Instead, the response of zooplankton to contacts with phytoplankton is likely to be delayed due to gestation. For example, in ^[35], the authors discussed Hopf bifurcation in the presence of time delay required for toxin-phytoplankton maturation. The universality of time-delay coupled system indicates its importance, applicability and practicability in a wide range of biological systems ^[36,37]. In fact, time delay may change the qualitative behavior of dynamic system ^[38,39].

In ^[40], the authors considered a fractional nutrient-phytoplankton-zooplankton system as follows

Based on the above model, we classify phytoplankton into non-toxic phytoplankton and toxic phytoplankton and put forward an improved fractional order four-dimensional ecological epidemiological model with delay. The system is established as follows:

subjected to the biologically feasible initial condition:

where $\phi(t)$, $\psi(t)$ and $\zeta(t)$ are continuous function defined on $t\in[-\tau, 0]$.

The meaning of state variables and parameters are listed in Table 1; $D^{\alpha}\ (0 < \alpha < 1)$ denotes Caputo fractional differential operator, and the model are based on the following scenarios:

(H1) $X(t)$, $Y_{1}(t)$, $Y_{2}(t)$ and $Z(t)$ represent nutrient population, phytoplankton population, toxic phytoplankton population and zooplankton population, respectively.

(H2) In real ecosystems, phytoplankton compete with each other for essential resources: nutrients and light. So as the model in the ^[41], we assume that, for nutrient $X(t)$, phytoplankton population $Y_{1}(t)$ is in competition with toxic phytoplankton population $Y_{2}(t)$, $h_{1}$ and $h_{2}$ represent the influence on $Y_{1}(t)$ and $Y_{2}(t)$ in the competition, respectively.

(H3) Zooplankton do not grow instantaneously after consuming phytoplankton, and pregnancy of predators requires a discrete time delay $\tau$.

(H4) Zooplankton populations feed only on phytoplankton, and only some of the dead phytoplankton and zooplankton are recycled into nutrients.

(H5) The toxic phytoplankton has both positive and negative effect on zooplankton, corresponding to the term $\theta_{2}\eta_{2}Y_{2}(t-\tau)Z(t-\tau)$ and $-\delta Y_{2}(t)Z(t)$ in the last equation of the system(1.2).

The above assumption (H5) is based on the result in ^[42]. In fact, the authors concluded that the food selection mechanism of plankton may not yet mature. When different algae of the same species (toxic and non-toxic) coexist, some zooplankton may have poor ability to select between toxic and non-toxic algae, and even show a slight preference for toxic strains.

The present paper is organized as follows. In section 2, some preliminaries are presented. In section 3, qualitative analysis of the system is performed. In section 4, some numerical examples and simulations are exploited to verify the theoretical results. In the last section, some conclusions and discussions are provided.

2.   Preliminaries

For convenience, we list some of the basic definitions and lemmas of the fractional calculus. In fractional-order calculus, there are many fractional-order integration and fractional-order differentiation that have been defined, for example, the Grunwald-Letnikov (GL) definition, the Riemann-Liouville (RL) definition and the Caputo definition. Since the initial condition is the same as the form of integral differential equation, we will adopt the definition of Caputo in this paper.

Definition 2.1. ^[34] The Riemann-Liouville fractional integral of order $\alpha > 0$ for a function $f:R^{+}\rightarrow R$ is defined by

Based on this definition of Riemann-Liouville fractional integral, the fractional-order derivative in Riemann-Liouville sense and Caputo sense are given.

Definition 2.2. ^[34] The Riemann-Liouville fractional derivative of order $\alpha > 0$ for a function $f: R^{+}\rightarrow R$ is defined by

where $k-1\leq\alpha < k, \ k\in N$ and $\Gamma(\cdot)$ is the Gamma function, $\Gamma(\alpha) = \int_{0}^{+\infty}t^{\alpha-1}e^{-t}\mathrm{d}t.$

In particular, when $0 < \alpha < 1$, we have

Definition 2.3. ^[34] The Caputo fractional derivative of order $\alpha > 0$ for a function $f:R^{+}\rightarrow R$ is defined by

where $k-1\leq\alpha < k, \ k\in N$ and $f^{(m)}(t)$ is the m-order derivative of f(t). In particular, when $0 < \alpha < 1$, we have

Definition 2.4. ^[34] The two-parameter Mittag-Leffler function is defined by

When $\beta = 1$, the two-parameter Mittag-Leffler function becomes to the one-parameter Mittag- Leffler function, i.e.

Theorem 2.5. ^[43] Consider the following commensurate fractional-order system:

with $0 < \alpha < 1$ and $x\in R^{n}$. The equilibrium points of the above system are calculated by solving the equation: $f(x) = 0$. These points are locally asymptotically stable if all eigenvalues $\lambda_{i}$ of the Jacobian matrix evaluated at the equilibrium points satisfy the inequality: ${|\mathrm{arg}(\lambda_{i})| > \frac{\alpha \pi}{2}}$.

3.   Qualitative analysis of the system (1.2)

3.1.   The existence of equilibriums

Since the proof of the positivity and boundedness of the solution of the system(1.2) is similar to Theorem 2 and Theorem 3 in the Ref.^[38], we will not prove it here.

The equilibriums of model (1.2) are obtained by solving the following algebraic system

By simple calculation, we obtain seven equilibriums of system (1.2), namely:

(1) $E_{0} = (X^{(0)}, \ 0, \ 0, \ 0)$ with $X^{(0)} = {\frac{\Lambda}{\mu}}$.

(2) $E_{1} = (X^{(1)}, \ Y^{(1)}_{1}, \ 0, \ 0)$ with $X^{(1)} = {\frac{\mu_{1}}{k_{1}b_{1}}}, \ Y^{(1)}_{1} = {\frac{\mu(1-R_{1})}{b_{1}(\frac{c_{1}R_{1}}{b_{1}X^{(0)}}-1)}}$, where ${R_{1} = \frac{X^{(0)}}{X^{(1)}}}$. And the feasibility conditions for $E_{1}$ are simplified as:

(3) $E_{2} = (X^{(2)}, \ 0, \ Y^{(2)}_{2}, \ 0)$ with $X^{(2)} = {\frac{\mu_{2}}{k_{2}b_{2}}}, \ Y^{(2)}_{2} = {\frac{\mu(1-R_{2})}{b_{2}(\frac{c_{2}R_{2}}{b_{2}X^{(0)}}-1)}}$, where ${R_{2} = \frac{X^{(0)}}{X^{(2)}}}$. And the feasibility conditions for $E_{2}$ are simplified as:

(4) $E_{3} = (X^{(3)}, \ Y_{1}^{(3)}, \ Y_{2}^{(3)}, \ 0)$ with $Y_{1}^{(3)} = {\frac{k_{2}b_{2}X^{(3)}-\mu_{2}}{h_{_{2}}}}$, $Y_{2}^{(3)} = {\frac{k_{1}b_{1}X^{(3)}-\mu_{1}}{h_{_{1}}}}$, and $X^{(3)}$ is uniquely determined by the following equation:

where

If ${\Lambda > \frac{\mu_{2}c_{1}}{h_{2}}+\frac{\mu_{1}c_{2}}{h_{1}}}$, then Descartes rule of sign ensures that the above Eq.(3.2) possesses a uniquely positive root. And the feasibility conditions for $E_{3}$ are simplified as:

(5) $E_{4} = (X^{(4)}, \ Y_{1}^{(4)}, \ 0, \ Z^{(4)})$ with

The feasibility conditions of $E_{4}$ are simplified as:

(6) $E_{5} = (X^{(5)}, \ 0, \ Y_{2}^{(5)}, \ Z^{(5)})$ with

Considering the biological background, we assume $\theta_{2}\eta_{2} > \delta$ is reasonable, and the feasibility conditions for $E_{5}$ are simplified as:

(7) $E_{6} = (X^{(6)}, \ Y_{1}^{(6)}, \ Y_{2}^{(6)}, \ Z^{(6)})$ with

$Z^{(6)} = {\frac{(\theta_{1}\eta_{1}h_{1}k_{2}b_{2}+\theta_{2}\eta_{2}h_{2}k_{1}b_{1}-\delta h_{2}k_{1}b_{1})X^{(6)}+\delta h_{2}\mu_{1}-\theta_{1}\eta_{1}h_{1}\mu_{2}-\theta_{2}\eta_{2}h_{2}\mu_{1}-\mu_{3}h_{1}h_{2}}{\eta_{1}(\theta_{1}\eta_{2}h_{1}+ \theta_{2}\eta_{2}h_{2}-\delta h_{2})}}$. And $X^{(6)}$ is uniquely determined by the following equation:

where

If $b_{1}b_{3} < 0$, then Descartes rule of sign ensures that the above Eq.(3.3) possesses a uniquely positive root. The feasibility conditions for $E_{6}$ are simplified as: $Y_{1}^{(6)}, Y_{2}^{(6)}, Z^{(6)} > 0, b_{1}b_{3} < 0.$

Remark 3.1. (1) The necessary conditions for the existence of $E_{3}$ are $R_{1} > 1$ and $R_{2} > 1$.

(2) Because of the complexity of computation, we have not obtain the exact formula of positive equilibrium $E_{6}$.

3.2.   Local stability of Equilibriums

In this subsection we discuss the stability of each equilibrium when $\tau$ = 0.

Obviously, the eigenvalues of the Jacobian matrix of system (1.2) at equilibrium $E_{0}$ are $\lambda_1 = -\mu < 0, \ \lambda_2 = -\mu_{3} < 0, \ \lambda_3 = k_{1}b_{1}X^{(0)}-\mu_{1}, \ \lambda_4 = k_{2}b_{2}X^{(0)}-\mu_{2},$ so we get the following result.

Theorem 3.2. If $R_{1} < 1$ and $R_{2} < 1$, then the disease-free equilibrium $E_{0}$ is locally asymptotically stable and it is unstable if $R_{1} > 1$ or $R_{2} > 1$.

The Jacobian matrix of system (1.2) at equilibrium $E_{1}$ is

The characteristic equation at the equilibrium $E_{1}$ is

Thus, we get the following result.

Theorem 3.3. If $1 < R_{1} < \frac{X^{(0)}b_{1}}{c_{1}}$, $R_{2} < R_{1}$ and $\mu_{3} > \theta_{1}\eta_{1}Y^{(1)}_{1}$, then the equilibrium $E_{1}$ is locally asymptotically stable.

The Jacobian matrix of system (1.2) at equilibrium $E_{2}$ is

The characteristic equation at the equilibrium $E_{2}$ is

So we get the following result.

Theorem 3.4. If $1 < R_{2} < \frac{X^{(0)}b_{2}}{c_{2}}$, $R_{2} > R_{1}$ and $\mu_{3} > (\theta_{2}\eta_{2}-f)Y^{(2)}_{2}$, then the equilibrium $E_{2}$ is locally asymptotically stable.

The Jacobian matrix of system (1.2) at equilibrium $E_{3}$ is

with the characteristic equation

where

Denote $D(P)$ denote the discriminant of a polynomial

$Q(\lambda) = \lambda^{3}+A_{1}\lambda^{2}+A_{2}\lambda+A_{3}$. Then

$D(Q) = 18A_{1}A_{2}A_{3}+(A_{1}A_{2})^{2}-4A_{3}A_{1}^{3}-4A_{2}^{3}-27A_{3}^{3}$.

If $\mu_{3} > \theta_{1}\eta_{1}Y_{1}^{(3)}+\theta_{2}\eta_{2}Y_{2}^{(3)}-\delta Y_{2}^{(3)}$, in order to discuss the stability of the equilibrium $E_{3}$, we get the following result by use of the same method as in Ref ^[44].

Proposition 3.5. The equilibrium $E_{3}$ is asymptotically stable if one of the following conditions holds for polynomial $Q$ and $D(Q)$:

$(1)D(Q) > 0, A_{1} > 0, A_{3} > 0$ and $A_{1}A_{2} > A_{3}$.

$(2)D(Q) < 0, A_{1}\geq0, A_{2}\geq0, A_{3}\geq0$ and $\alpha < \frac{2}{3}$.

The Jacobian matrix of system (1.2) at equilibrium $E_{4}$ is

where

The characteristic equation at the equilibrium $E_{4}$ is

where

Denote $D(P)$ denote the discriminant of a polynomial

$Q(\lambda) = \lambda^{3}+B_{1}\lambda^{2}+B_{2}\lambda+B_{3}$. Then

$D(Q) = 18B_{1}B_{2}B_{3}+(B_{1}B_{2})^{2}-4B_{3}B_{1}^{3}-4B_{2}^{3}-27B_{3}^{3}$.

If $\mu_{2} > k_{2}b_{2}X^{(4)}-\eta_{2}Z^{(4)}-h_{2}Y_{1}^{(4)}$, in order to discuss the stability of the equilibrium $E_{4}$, we get the following result by use of the same method as in Ref ^[44].

Proposition 3.6. The equilibrium $E_{4}$ is asymptotically stable if one of the following conditions holds for polynomial $Q$ and $D(Q)$:

$(1)D(Q) > 0, B_{1} > 0, B_{3} > 0$ and $B_{1}B_{2} > B_{3}$.

$(2)D(Q) < 0, B_{1}\geq0, B_{2}\geq0, B_{3}\geq0$ and $\alpha < \frac{2}{3}$.

The Jacobian matrix of system (1.2) at equilibrium $E_{5}$ is

where

The characteristic equation at the equilibrium $E_{5}$ is

where

Denote $D(P)$ denote the discriminant of a polynomial

$Q(\lambda) = \lambda^{3}+C_{1}\lambda^{2}+C_{2}\lambda+C_{3}$. Then

$D(Q) = 18C_{1}C_{2}C_{3}+(C_{1}C_{2})^{2}-4C_{3}C_{1}^{3}-4C_{2}^{3}-27C_{3}^{3}$.

If $\mu_{1} > k_{1}b_{1}X^{(5)}-\eta_{1}Z^{(5)}-h_{1}Y_{2}^{(5)}$, in order to discuss the stability of the equilibrium $E_{5}$, we get the following result by use of the same method as in Ref ^[44].

Proposition 3.7. The equilibrium $E_{5}$ is asymptotically stable if one of the following conditions holds for polynomial $Q$ and $D(Q)$:

$(1)D(Q) > 0, C_{1} > 0, C_{3} > 0$ and $C_{1}C_{2} > C_{3}$.

$(2)D(Q) < 0, C_{1}\geq0, C_{2}\geq0, C_{3}\geq0$ and $\alpha < \frac{2}{3}$.

Theorem 3.8. The equilibrium $E_{6}$ is locally asymptotically stable if the following conditions hold:

(H5) $X^{(6)} > \max {\left\{ \frac{c_{1}}{b_{1}}, \ \frac{c_{2}}{b_{2}}\right\}}$,

(H6) ${\frac{k_{1}b_{1}}{\eta_{1}} < \frac{k_{2}b_{2}}{\eta_{2}}}$,

(H7) ${k_{3} > \frac{\eta_{2}(\mu+b_{1}Y_{1}^{(6)}+b_{2}Y_{2}^{(6)})}{c_{3}b_{2}}}$,

(H8) ${(e_{1}e_{2}-e_{3})e_{3}-e_{4}e_{1}^{2} > 0}$.

Proof. The Jacobian matrix of system (1.2) at equilibrium $E_{6}$ is

where

The characteristic equation at the equilibrium $E_{6}$ is

where

By simple calculation, if $X^{(6)} > \max {\left\{ \frac{c_{1}}{b_{1}}, \ \frac{c_{2}}{b_{2}}\right\}}$, then $e_{1}e_{2} > e_{3}$; if ${k_{3} > \frac{\eta_{2}(\mu+b_{1}Y_{1}^{(6)}+b_{2}Y_{2}^{(6)})}{c_{3}b_{2}}}$ and ${\frac{k_{1}b_{1}}{\eta_{1}} < \frac{k_{2}b_{2}}{\eta_{2}}}$, then $e_{4} > 0$. In summary, the condition of the Routh-Hurwitz criterion above is satisfied for Eq.(3.9), that is,

hold. So, all the roots of this Eq.(3.9) have negative real part. This ends the proof.

Remark 3.9 In Theorem (3.8), (H5)-(H8) is a sufficient condition for equilibrium $E_{6}$ to be stable, and the necessary and sufficient condition for $E_{6}$ to be stable is all roots of Eq.(3.9) satisfy $|\mathrm{arg}(\lambda_{i})| > \frac{\alpha\pi}{2}$.

3.3.   Hopf bifurcation

In this subsection, according to the research methods in literature ^[22,38,45], we study the Hopf bifurcation with time delay as the parameter.

we will analyze the Hopf bifurcation of $E_{6}$ when $\tau > 0$, and the characteristic equation at the equilibrium $E_{6}$ is

where

Assume that $s = i\omega = \omega(\cos\frac{\pi}{2}+i\sin\frac{\pi}{2})$, $\omega > 0$ is a root of Eq.(3.10).

Substituting $s = i\omega$ into Eq.(3.10), one gets

and separating the real and imaginary parts of it, it results in

$R_{i}, I_{i}$ are defined as follows:

It can be acquired from Eq. (3.12) that

Adding the squares of the two equations of Eq.(3.12), we obtain

where $M$ is a polynomial containing $\omega^{7\alpha}, \ \omega^{6\alpha}, \ \omega^{5\alpha}, \ \omega^{4\alpha}, \ \omega^{3\alpha}, \ \omega^{2\alpha}, \ \omega^{\alpha}$, and $N$ is a constant.

Define

Suppose that $N < 0$. Thus, $h(\omega)$ has at least one positive root. The delay $\tau$ is regarded as a bifurcation parameter. Let $s(\omega) = \xi(\tau)+i\omega(\tau)$ be the Eq.(3.10) such that for some initial value of the bifurcation parameter $\tau_{0}$ we have $\xi(\tau_{0}) = 0$, $\omega(\tau_{0}) = \omega_{0}$. Without loss of generality, we assume $\omega(0) > 0$. From Eq.(3.13), one can conclude

where

To derive the condition of the occurrence for Hopf bifurcation, we have the following Lemma.

Lemma 3.10. Assume that $N < 0$, then Hopf bifurcation occurs provided $h^{'}(\omega_{0})\neq 0$.

Proof. Differentiating both sides of Eq.(3.10) with respect to $\tau$, it can be obtained that

Hence, one gets

Substitute $s = i\omega_{0}$ into Eq.(3.17), we have

where

then $\mathrm{sign}\left\{\frac{\mathrm{d}\mathrm{Re}(\lambda)}{\mathrm{d}\tau}|_{\tau = \tau_{0}}\right\} = \mathrm{sign}\left\{\mathrm{Re}\left[\left(\frac{\mathrm{d}\lambda} {\mathrm{d}\tau}\right)^{-1}|_{\tau = \tau_{0}}\right]\right\} = \mathrm{sign}\{h^{'}(\omega_{0})\}$.

Obviously, if $h^{'}(\omega_{0})\neq 0$ the transversality condition holds, and Hopf bifurcation occurs at $\tau = \tau_{0}$.

Theorem 3.11. Suppose that (H5)-(H8) and $N < 0$ hold, then the positive equilibrium $E_{6}$ of system (1.2) is asymptotically stable when $\tau\in[0, \tau_{0})$, $h^{'}(\omega_{0}) < 0$ and unstable when $\tau > \tau_{0}$, $h^{'}(\omega_{0}) > 0$. When $\tau = \tau_{0}$, $h^{'}(\omega_{0})\neq0$ a Hopf bifurcation occurs, that is a family of periodic solutions bifurcates from $E_{6}$ as $\tau$ passes through the critical value $\tau_{0}$.

4.   Numerical simulations

In this section, some numerical examples are presented to verify the theoretical results. The simulation are based on Adama-Bashforth-Moulton predictor-corrector scheme ^[46].

Example 1. For the following set of parameters: $\Lambda =$ 1, $b_{1} =$ 0.3, $b_{2} =$ 0.25, $c_{1} =$ 0.06, $c_{2} =$ 0.06, $c_{3} =$ 0.06, $k_{1} =$ 0.7, $k_{2} =$ 0.7, $\eta_{1} =$ 2.1, $\eta_{2} =$ 0.2, $\mu =$ 1, $\mu_{1} =$ 0.5, $\mu_{2} =$ 0.3, $\mu_{3} =$ 0.3, $\theta_{1} =$ 0.6, $\theta_{2} =$ 0.7, $\delta =$ 0.1, $h_{1} =$ 0.2, $h_{2} =$ 0.1.

In this case $R_{1} = 0.42 < 1$, $R_{2} = 0.5833 < 1$. From Figure 1, we can see that the equilibrium $E_{0} = (1, \ 0, \ 0, \ 0)$ is stable for different values of $\alpha$ and different sets of initial values: $\left[X(0), \ Y_{1}(0), \ Y_{2}(0), \ Z(0)\right] = \left[0.9, \ 0.2, \ 0.2, \ 0.2\right]$, $\left[1.2, \ 0.5, \ 0.3, \ 0.6 \right]$.

Example 2. For the following set of parameters: $\Lambda =$ 0.5, $b_{1} =$ 2.5, $b_{2} =$ 0.3, $c_{1} =$ 0.01, $c_{2} =$ 0.01, $c_{3} =$ 0.06, $k_{1} =$ 0.5, $k_{2} =$ 0.7, $\eta_{1} =$ 0.2, $\eta_{2} =$ 2.1, $\mu =$ 1, $\mu_{1} =$ 0.3, $\mu_{2} =$ 0.5, $\mu_{3} =$ 0.3, $\theta_{1} =$ 0.5, $\theta_{2} =$ 0.1, $\delta =$ 0.01, $h_{1} =$ 0.2, $h_{2} =$ 0.1.

In this case, $1 < R_{1} = 2.0833 < \frac{X^{(0)}b_{1}}{c_{1}}$, $R_{2} = 0.21 < R_{1}$ and $\mu_{3} > \theta_{1}\eta_{1}Y^{(1)}_{1}$. From Figure 2, we can see that the equilibrium $E_{1} = (0.24, \ 0.4407, \ 0, \ 0)$ is stable for different values of $\alpha$ and different sets of initial values: $[X(0), \ Y_{1}(0), \ Y_{2}(0), \ Z(0)] = [0.1, 0.6, \ 0.5, \ 0.2]$, $[0.2, 0.5, \ 0.8, \ 0.6]$.

Example 3. For the following set of parameters: $\Lambda =$ 0.5, $b_{1} =$ 0.3, $b_{2} =$ 2.5, $c_{1} =$ 0.01, $c_{2} =$ 0.01, $c_{3} =$ 0.06, $k_{1} =$ 0.5, $k_{2} =$ 0.6, $\eta_{1} =$ 2.1, $\eta_{2} =$ 0.2, $\mu =$ 1, $\mu_{1} =$ 0.3, $\mu_{2} =$ 0.5, $\mu_{3} =$ 0.5, $\theta_{1} =$ 0.1, $\theta_{2} =$ 0.1, $\delta =$ 0.01, $h_{1} =$ 0.2, $h_{2} =$ 0.1.

In this case, $1 < R_{2} = 1.5 < \frac{X^{(0)}b_{2}}{c_{2}}$, $R_{2} > R_{1}$ and $\mu_{3} > (\theta_{2}\eta_{2}-f)Y^{(2)}_{2}$. From Figure 3, we can see that the equilibrium $E_{2} = (0.3333, \ 0, \ 0.2024, \ 0)$ is stable for different values of $\alpha$ and different sets of initial values: $[X(0), \ Y_{1}(0), \ Y_{2}(0), \ Z(0)] = [0.2, \ 0.1, \ 0.5, \ 0.2]$, $[0.1, \ 0.05, \ 0.3, \ 0.6]$.

Remark 4.1. The above three examples corresponding to the following ecological interpretation.

(1) Figure 1 indicates that if $R_{1} < 1$, $R_{2} < 1$, then the phytoplankton can not survive and the zooplankton will also die out. However, this phenomenon usually does not happen in the real world.

(2) Figure 2 indicates that if $R_{1} > 1$, $R_{1} > R_{2}$, then the non-toxic phytoplankton will win the competition between phytoplankton and toxic phytoplankton, while the zooplankton will die out due to excessive mortality.

(3) Figure 3 indicates that if $R_{2} > 1$, $R_{2} > R_{1}$, then the toxic phytoplankton will win the competition between phytoplankton and toxic phytoplankton, while the zooplankton will die out due to excessive mortality.

Example 4. For the following set of parameters: $\Lambda =$ 1.4, $b_{1} =$ 2.8, $b_{2} =$ 2.3, $c_{1} =$ 0.01, $c_{2} =$ 0.01, $c_{3} =$ 0.06, $k_{1} =$ 0.95, $k_{2} =$ 0.1, $\eta_{1} =$ 2.1, $\eta_{2} =$ 0.26, $\mu =$ 0.2, $\mu_{1} =$ 0.2, $\mu_{2} =$ 0.8, $\mu_{3} =$ 0.5, $\theta_{1} =$ 0.2, $\theta_{2} =$ 0.5, $\delta =$ 0.01, $h_{1} =$ 0.1, $h_{2} =$ 0.1.

In this case, simple calculation indicates that the sufficient condition (2) of Proposition 3.6 is satisfied. From Figure 4 we can see that the equilibrium $E_{4} = (0.4067, \ 1.1905, \ 0, \ 0.4199)$ is stable for different values of $\alpha$ and different sets of initial values: $[X(0), \ Y_{1}(0), \ Y_{2}(0), \ Z(0)] = [0.2, \ 0.3, \ 0.5, \ 0.2]$, $[0.3, \ 0.4, \ 0.8, \ 0.6]$.

Example 5. For the following set of parameters: $\Lambda =$ 1.4, $b_{1} =$ 2.3, $b_{2} =$ 2.5, $c_{1} =$ 0.01, $c_{2} =$ 0.01, $c_{3} =$ 0.06, $k_{1} =$ 0.1, $k_{2} =$ 0.95, $\eta_{1} =$ 0.26, $\eta_{2} =$ 2.1, $\mu =$ 0.2, $\mu_{1} =$ 0.8, $\mu_{2} =$ 0.2, $\mu_{3} =$ 0.5, $\theta_{1} =$ 0.5, $\theta_{2} =$ 0.2, $\delta =$ 0.01, $h_{1} =$ 0.1, $h_{2} =$ 0.1.

In this case, simple calculation indicates that the sufficient condition (2) of Proposition 3.7 is satisfied. From Figure 5 we can see that the equilibrium $E_{5} = (0.4422, \ 0, \ 1.2195, \ 0.4048)$ is stable for different values of $\alpha$ and different sets of initial values: $[X(0), \ Y_{1}(0), \ Y_{2}(0), \ Z(0)] = [0.2, \ 0.5, \ 0.3, \ 0.2]$, $[0.3, \ 0.8, \ 0.4, \ 0.6]$.

Remark 4.2. The above examples corresponding to the following ecological interpretation.

(1) Figure 4 indicates that if the toxic phytoplankton is less competitive than the non-toxic phytoplankton, then the toxic phytoplankton will die out.

(2) Figure 5 indicates that if the toxic phytoplankton win the competition between non-toxic phytoplankton and toxic phytoplankton, then the zooplankton may still survive under certain conditions, that is, nutrients, toxic phytoplankton and zooplankton may theoretically coexist. However, this phenomenon usually does not appear in real world.

(3) From Figures 1–5 we can see that as the value of $\alpha$ decreases, the steady speed becomes slow for each equilibrium. This indicates that the value of $\alpha$ has obvious effects on the dynamical behaviors of the system.

Example 6. For the following set of parameters: $\alpha =$ 0.8, $\Lambda =$ 1.4, $b_{1} =$ 0.32, $b_{2} =$ 0.54, $c_{1} =$ 0.06, $c_{2} =$ 0.08, $c_{3} =$ 0.6, $k_{1} =$ 0.7, $k_{2} =$ 0.6, $\eta_{1} =$ 1.8, $\eta_{2} =$ 0.6, $\mu =$ 0.2, $\mu_{1} =$ 0.4, $\mu_{2} =$ 0.9, $\mu_{3} =$ 0.5, $\theta_{2} =$ 0.5, $\delta =$ 0.1, $h_{1} =$ 0.1, $h_{2} =$ 0.1.

In this example, we will consider the influence of toxic, i.e., $\theta_{1}$. Here, we choose $\theta_{1} = 0, 0.6$, with the initial conditions: $[X(0), \ Y_{1}(0), \ Y_{2}(0), \ Z(0)]$ = $[0.2, \ 0.2, \ 0.3, \ 0.2], \ [0.8, \ 0.8, \ 0.6, \ 0.8]$.

Example 7. In this example, we will consider the influence of $\alpha$.

(1) For the following set of parameters: $\Lambda =$ 1.4, $b_{1} =$ 2.3, $b_{2} =$ 2.5, $c_{1} =$ 0.01, $c_{2} =$ 0.01, $c_{3} =$ 0.06, $k_{1} =$ 0.1, $k_{2} =$ 0.95, $\eta_{1} =$ 0.26, $\eta_{2} =$ 2.1, $\mu =$ 0.2, $\mu_{1} =$ 0.8, $\mu_{2} =$ 0.2, $\mu_{3} =$ 0.5, $\theta_{1} =$ 0.5, $\theta_{2} =$ 0.2, $\delta =$ 0.01, $h_{1} =$ 0.1, $h_{2} =$ 0.1, with the initial conditions: $[X(0), \ Y_{1}(0), \ Y_{2}(0), \ Z(0)]$ = $[0.2, \ 0.5, \ 0.3, \ 0.2]$.

(2) For the following set of parameters: $\Lambda =$ 1.4, $b_{1} =$ 0.32, $b_{2} =$ 0.54, $c_{1} =$ 0.06, $c_{2} =$ 0.08, $c_{3} =$ 0.6, $k_{1} =$ 0.7, $k_{2} =$ 0.6, $\eta_{1} =$ 1.8, $\eta_{2} =$ 0.6, $\mu =$ 0.2, $\mu_{1} =$ 0.4, $\mu_{2} =$ 0.9, $\mu_{3} =$ 0.5, $\theta_{1} =$ 0.6, $\theta_{2} =$ 0.5, $\delta =$ 0.1, $h_{1} =$ 0.1, $h_{2} =$ 0.1, with the initial conditions: $[X(0), \ Y_{1}(0), \ Y_{2}(0), \ Z(0)]$ = $[0.5, \ 0.5, \ 0.5, \ 0.2]$.

Remark 4.3. (1) From Figure 7, we can see that if the value of $\alpha$ is relatively big(i.e., $\alpha$ = 1, 0.8), then the equilibrium $E_5$ is locally stable; if the value of $\alpha$ is relatively small(i.e., $\alpha$ = 0.4), then the equilibrium is unstable, and oscillation may occur.

(2) From Figure 8, we can see that if the value of $\alpha$ is relatively big(i.e., $\alpha$ = 1, 0.7), then the equilibrium $E_6$ is locally stable; if the value of $\alpha$ is relatively small(i.e., $\alpha$ = 0.1), then the equilibrium $E_0$ is locally stable.

(3) Figure 7 and Figure 8 indicate that if the value of $\alpha$ is relatively small, the system will be destabilized.

Example 8. For the following set of parameters: $\alpha =$ 0.89, $\Lambda =$ 1.4, $b_{1} =$ 0.32, $b_{2} =$ 0.54, $c_{1} =$ 0.06, $c_{2} =$ 0.08, $c_{3} =$ 0.6, $k_{1} =$ 0.7, $k_{2} =$ 0.6, $\eta_{1} =$ 1.8, $\eta_{2} =$ 0.6, $\mu =$ 0.2, $\mu_{1} =$ 0.4, $\mu_{2} =$ 0.9, $\mu_{3} =$ 0.5, $\theta_{1} =$ 0.6, $\theta_{2} =$ 0.5, $\delta =$ 0.1, $h_{1} =$ 0.1, $h_{2} =$ 0.1.

In this example, the equilibrium $E_{6} = (3.2026, \ 0.4114, \ 0.2784, \ 0.1608)$, and our main aim is to study the effect of time delay on the stability of the system.

(1) From Figure 9, we can see that if $\tau = 0$, then the equilibrium $E_{6}$ is stable with different sets of initial value: $[X(0), \ Y_{1}(0), \ Y_{2}(0), \ Z(0)] = [0.5, \ 0.5, \ 0.5, \ 0.2]$, $[2.8, \ 1.8, \ 3.3, \ 1.3]$.

(2) From Figure 10, we will see that if $\tau = 3 < \tau_{0}\approx4.4671$, then the equilibrium $E_{6}$ is stable with different sets of initial value: $[X(0), \ Y_{1}(0), \ Y_{2}(0), \ Z(0)] = [0.5, \ 0.5, \ 0.5, \ 0.2]$, $[2.8, \ 1.8, \ 3.3, \ 1.3]$.

(3) From Figure 11, we can see that if $\tau = 4.2 < \tau_{0}\approx4.4671$, then the equilibrium $E_{6}$ is stable with different sets of initial value: $[X(0), \ Y_{1}(0), \ Y_{2}(0), \ Z(0)] = [0.5, \ 0.5, \ 0.5, \ 0.2]$, $[2.8, \ 1.8, \ 3.3, \ 1.3]$.

(4) From Figure 12, we will see that if $\tau = 5 > \tau_{0}\approx4.4671$, then periodic oscillation occurs, and the equilibrium $E_{6}$ will lose its stability and periodic solutions appear through Hopf bifurcation, with different sets of initial value: $[X(0), \ Y_{1}(0), \ Y_{2}(0), \ Z(0)] = [0.5, \ 0.5, \ 0.5, \ 0.2]$, $[0.8, \ 0.8, \ 0.3, \ 0.3]$.

(5) From Figure 13, we can see that if $\tau = 6 > \tau_{0}\approx4.4671$, then similar periodic oscillation occurs as that in Figure 12, but with larger amplitude.

Remark 4.4. From the above example, we can see that the time delay has an effect of destabilizing the equilibrium $E_{6}$. In other words, the larger the value of time delay is, the more possible that the equilibrium $E_{6}$ lose its stability.

5.   Discussion

In this paper, a fractional-order mathematical model is constructed to describe the active of nutrient-phytoplankton-toxic phytoplankton-zooplankton.

Through qualitative analysis, we get the following results.

$\lozenge$ We figure out the sufficient conditions for the existence and local stability of $E_{0}$, $E_{1}$, $E_{2}$, $E_{3}$, $E_{4}$, $E_{5}$, $E_{6}$ for $\tau = 0$.

$\lozenge$ By using time delay as a bifurcation parameter, the existence of Hopf bifurcation is analyzed in detail. We find that if $\tau < \tau_{0}$, then the equilibrium $E_{6}$ is locally stable; while it is unstable if $\tau > \tau_{0}$ and Hopf bifurcation may occur near $\tau_{0}$.

Through numerical simulation we get the following results.

$\lozenge$ Figure 1 shows the stability of equilibrium $E_{0}$ for different values of $\alpha$.

$\lozenge$ Figure 2 shows the stability of equilibrium $E_{1}$ for different values of $\alpha$.

$\lozenge$ Figure 3 shows the stability of equilibrium $E_{2}$ for different values of $\alpha$.

$\lozenge$ Figure 4 shows the stability of equilibrium $E_{4}$ for different values of $\alpha$.

$\lozenge$ Figure 5 shows the stability of equilibrium $E_{5}$ for different values of $\alpha$.

$\lozenge$ Figure 6 shows the effect of parameter $\theta_{1}$ on the system (1.2).

$\lozenge$ Figure 7 and Figure 8 indicate that the value of $\alpha$ is closely relate to the stability of each equilibrium. The stability of each equilibrium becomes weaker as the value of $\alpha$ decreases.

$\lozenge$ Figures 9–11 show that $E_{6}$ is stable if $\tau\in[0, \tau_{0})$; Figures 12–13 show that $E_{6}$ is unstable if $\tau > \tau_{0}$ and periodic oscillation may occur; Figures 9–13 indicate that the impact of $\tau$ on the dynamics of the system is crucial.

$\lozenge$ Table 2 and Figure 14 show that the value of $\tau_{0}$ arises as the value of $\alpha$ increases.

Remark 5.1. When $\tau = 0$, $Y_{1} = Y_{2}$ and $h_{1} = h_{2} = 0$, then system (1.2) will degenerated to the model in ^[40].

In the model of this paper, phytoplankton is divided into two class, non-toxic phytoplankton and toxic phytoplankton. Some experimental data suggested that some zooplankton are capable of selecting for nontoxic phytoplankton, a mechanism that allows toxic phytoplankton to coexist with nontoxic phytoplankton ^[47]. This is also consistent with the results in Figure 9. Although some zooplankton have the ability to distinguish between toxic and non-toxic plants, some other experiments have shown that some zooplankton might not be able to distinguish between toxic and nontoxic algae, even shows a slight preference for toxic strains ^[48,49]. We can find from Figure 6 that once non-toxic phytoplankton become scarce, the zooplankton start to eat the toxic phytoplankton, even if the toxicity is weak, the zooplankton may become extinct. This is dangerous for the ecosystem.

Remark 5.2. Any ecosystem depends on the natural environment, and in the real natural environment there are more or less physical and chemical signals that interact with the ecosystem. This motivate us to consider stochastic effects to the ecosystem, that is, white noise should be included into the system. We leave this as our next work.

6.   Authors contributions

Each of the authors, Ruiqing Shi, Jianing Ren, Cuihong Wang contributed to each part of this work equally and read and approved the final version of the manuscript.

7.   Acknowledgements

This work is partly supported by National Natural Science Foundation of China (No. 61907027). The authors would like to thank the anonymous reviewers for their helpful comments, which improved the quality of this paper greatly.

Conflict of interest

The authors declare that they have no financial or non-financial competing interests.