Global analysis and optimal harvesting for a hybrid stochastic phytoplankton-zooplankton-fish model with distributed delays

1.   Introduction

In aquatic ecosystems, phytoplankton are taken as basic food source and the first trophic level while zooplankton are primary consumers of phytoplankton in food chains ^[1]. Some plankton models composed of phytoplankton and zooplankton were formulated, and dynamical behaviors of those models were investigated in the past two decades (e.g., ^{[2,3,4,5,6,7,8]}). Besides plankton, fish is an essential part in aquatic environments like fishponds, lakes, rivers, oceans, etc.. According to an experiment given in ^[9], it was clearly showed that the addition of fish to the chain of phytoplankton-zooplankton caused the reduction in the algae intake of zooplankton and the rapid growth of phytoplankton. Thus, it is meaningful to incorporate plankton-feeding fish into the plankton model to form the food chain relationship with three species involving fish, phytoplankton and zooplankton. Actually, the study of asymptotical behavior for the plankton-fish model is closely related to the sustainable development of aquatic ecosystems. Scheffe ^[10] originally accounted for the effects of planktivorous fish in the phytoplankton-zooplankton interaction model, and Malchow et al.^[11] have extended the model to a spatial one. Recently, Prabir Panja et al.^[12] formulated a toxin-producing phytoplankton-zooplankton-fish model and obtained some sufficient conditions on stability, the existence of equilibrium and bifurcation of the model. Amit Sharma et al. ^[13] proposed a delayed plankton-fish model with harvesting, and the bifurcation analysis of the system was carried out by taking the rate of harvesting as the bifurcation parameter. Meng and Wu ^[14] proposed a delayed phytoplankton-zooplankton-fish model with taxation and nonlinear fish harvesting, and gave the Hopf-bifurcation analysis for the model. Wei et al.^[15] considered a time-varying phytoplankton-zooplankton-fish system, and gave some sufficient conditions ensuring global asymptotical stability. Since many omnivorous fishes (like crucian carp, grass carp, carp, engraulis japonicus, cockerel, etc.) feed on both phytoplankton and zooplankton, there are more enriching and complex behaviors for the phytoplankton-zooplankton-fish system with the food chain relationship given in Figure 1 (e.g., ^[12,13]).

On the other hand, delay effects and environmental disturbances are unavoidable in the real world. Time delays occur frequently in many predator-prey models because predators can increase their quantity through digestion, absorption, reproduction and other processes after ingesting the food ^[16,17]. Meanwhile, environmental noises often affect dynamic behaviors of the population system ^[18,19,20] since the birth rate, death rate, environmental carrying rate and other system parameters of species are easily disturbed by noises. These noises usually involve white noise and telegraph noise. Telegraph noise can be regarded as switching without memory between two or more states, and the time spent in switching between two states is exponentially distributed ^[21,22]. Therefore, it is more realistic to formulate a phytoplankton-zooplankton-fish model with time delays and hybrid stochastic noises involving Brownian motion and Markov chain, which essentially belongs to the stochastic predator-prey system. In past decades, some interesting results on stochastic predator-prey models have been investigated. In ^[23], the authors studied the stationary distribution and global asymptotic stability of a three-species stochastic food-chain system without time-delay. In ^[24], the authors considered a three-species food chain stochastic system with a hidden Markov chain and proposed two kinds of special dissipative control strategies. One can refer to recent publications on dynamical behaviors of stochastic predator-prey systems ^{[25,26,27,28]}. To the best of our knowledge, however, there are few works to discuss the dynamics for the stochastic delayed phytoplankton-zooplankton-fish model with the food chain relationship given in Figure 1.

What's more, the optimal harvesting strategy is of great significance to the development of the ecosystem. In the phytoplankton-zooplankton-fish system, the eutrophication of the water body will be controlled by harvesting plankton ^[29]. Ones can directly benefit from fish, but overfishing may break the balance of the ecosystem. In response to the issue of resource sustainability, Clark and Mesterton-Gibbons et al. have established several types of predator-prey ecological models with optimal harvesting strategies, and discussed how to implement harvesting strategies to maintain fisheries sustainable development ^{[30,31,32,33]}. The optimal harvesting problem of the stochastic predator-prey model with time delays was investigated by the ergodic method in ^[34,35,36]. In the obtained optimization strategies, authors focused on the maximum of total species that has been harvested rather than the maximum of total economic income. By using Pontryagin's maximum principle, the optimal harvesting policy with the maximized present value of revenues was given for a deterministic phytoplankton-zooplankton model ^[3]. Nevertheless, there are few publications to investigate the optimal harvesting problem in economic income for the stochastic phytoplankton-zooplankton-fish model with time delays.

Motivated by the above discussion, we propose the following hybrid stochastic phytoplankton-zooplankton-fish system with distributed delays and harvesting

where $x_{i}(t), i = 1, 2, 3,$ is the population size of phytoplankton, zooplankton and fish at time $t$, respectively, $a_{1}(\cdot) > 0$ means the growth rate of species $x_{1}, a_{j}(\cdot) > 0 (j = 2, 3)$ stands for the death rate of species $x_{j}$, $c_{ii} > 0$ shows the intra-specific competition of the $i$th species, $c_{12}, c_{13}$ and $c_{23} > 0$ denote the capture rates, $c_{21}, c_{31}$ and $c_{32} > 0$ represent the conversation rate of the food, $h_{i}\ge 0$ is harvesting effort of phytoplankton, zooplankton and fish, respectively, $\tau_{ij}$ is time delay and $\mu_{ij}(\theta)$ is a deterministic and nondecreasing function defined on $[-\tau_{ij}, 0]$ satisfying $\int^{0}_{-\tau_{ij}}\mathrm{d}\mu_{ij}(\theta) = 1$, ${\beta_{i}(t)}_{t\ge 0}$ is standard independent Brownian motion defined on a complete probability space $\left(\Omega, \left\{\mathcal{F}_{t}\right\}_{t \geq 0}, P\right)$ with a filtration $\left\{\mathcal{F}_{t}\right\}_{t \geq 0}$ satisfying the usual conditions and $\sigma^2_i(\cdot)$ represents the intensity of the stochastic noise, $\{\gamma(t), t\ge0\}$ is a continuous-time Markov chain in a finite state space $\mathbb{S} = \{1, \cdots n\}$ with the generator $Q = (q_{ij})_{n\times n}$, which follows

and $\delta = 0$, $q_{ij}$ is the transition rate from $i$ to $j$ satisfying $q_{ij} > 0, q_{ii} = -\sum\limits_{j \ne i}q_{ij}.$

In this paper, we mainly study dynamical behaviors and optimal harvesting policy for the above phytoplankton-zooplankton-fish model with distributed delays, hybrid stochastic noises and harvesting. The main purpose and contribution of this paper are listed as follows. Firstly, we formulate a stochastic phytoplankton-zooplankton-fish model with distributed delays and harvesting, in which hybrid stochastic noises involve Brownian motion and Markov chain. Secondly, we give a global analysis of dynamics on persistence, extinction, stability and attractivity in terms of some system parameters for the stochastic delayed system. Lastly, we provide the optimal harvesting policy by solving the following optimization problem with the maximum of total economic income

where the harvesting effort $H$ is the decision variable, and the total profit $\Phi$ is the objective function, the unit profit $r_j = p_j-q_j > 0$, $p_j$ represents the unit market price of the species $x_j$, and $q_j$ is the cost unit price of harvesting the species $x_j, j = 1, 2, 3$, $W$ stands for the total fixed cost for harvesting three species.

The remaining part of this paper is organized as follows. In section 2, we give some definitions, assumptions and basic lemmas. In section 3, we give the global analysis of dynamic behavior on stability, persistence and extinction for the system (1). We obtain the sufficient and necessary condition for the optimal harvesting strategy and the maximum of harvesting yield in section 4. Lastly, we illustrate our main results in some examples and their simulations in section 5.

2.   Preliminary

In the model (1), we always suppose that Brownian motions $\beta_i(t)$ and the Markov chain $\gamma(t)$ are independent, and the Markov chain is irreducible. According to ^[37], $\gamma(t)$ is ergodic and has a unique stationary distribution $\xi = (\xi_{1}, \cdots \xi_{n})$ satisfying

For simplicity, we define the following notions.

$|C|$ represents the determinant of $C$, $C_{ij}$ stands for the complement minor of $|C|, i, j = 1, 2, 3$, for $m, n \in R$, $m \wedge n = \min\{m, n\}$, and let $I$ be the unit matrix.

Next, we will give the following assumptions.

${\bf (H_{1})}$ $|C| > 0, C_{i} > 0$, $\Delta_{j} > 0, j = 2, 3$;

${\bf (H_{2})}$ $\tilde C_{i} > 0, i = 1, 2, 3, C_{23} > 0, C_{12} > 0, C_{31} < 0$.

According to the above assumptions, we can see that the system (1) has a positive equilibrium state if there is no stochastic noise and harvesting and that species 1 and species $j, j = 2, 3$ can coexist without stochastic noise, harvesting and other predators ^[38].

For delayed stochastic system (1), the initial conditions are given in the following form

Here, we denote $\mathbb{C}\left([-\tau, 0]; \mathbb{R}_{+}^{3}\right)$ represents the family of all continuous functions from $[-\tau, 0] \rightarrow \left\{x \in \mathbb{R}^{3} | x_{1} > 0, x_{2} > 0, x_3 > 0 \right\}$ and $\tau = \max\limits_{i\neq j}\{\tau_{ij}\}$.

In the following, we give some definitions and lemmas to obtain our main results.

Definition 1. (see ^[39]) Let $X(t) = \left(x(t), \gamma(t)\right)\in \mathbb R^3_{+} \times \mathbb{S}$ be the solution of the system (1). Then for $i = 1, 2, 3,$

(a) the population $x_{i}(t)$ is said to be extinct if $\lim\limits_{t \to \infty}x_{i}(t) = 0$, a.s.;

(b) the population $x_{i}(t)$ is said to be stable in the mean if $\lim\limits_{t \to \infty}\frac{1}{t}\int_{0}^{t}x_{i}(s)\mathrm{d}s = c > 0$, a.s..

Definition 2. (see ^[40]) The system ({1) is said to be asymptotically stable in distribution if there exists a probability measure $\pi(\cdot \times \cdot)$ on $\mathbb{R}_{+}^{3} \times \mathbb{S}$ such that transition probability $p\left(t, \phi, i, \mathrm{d}y \times {j}\right)$ of the stochastic process $y(t)$ converges weakly to $\pi(dy \times {j})$ as $t \to \infty$ for every initial value $\phi \in C([-\tau, 0], {\mathbb R}^3_{+})$ with $\gamma(0) = i \in \mathbb{S}$.

For the biological significance, we first give the well-posedness and boundedness of the system (1).

Lemma 1. For any given initial value $\phi\in \mathbb{C}\left([-\tau, 0], \mathbb {R}_{+}^{3}\right)$ with $\gamma(0) = k \in \mathbb{S}$, the system (1) has a unique global positive solution $X(t) = (x(t), \gamma(t))\in \mathbb R^3_{+} \times \mathbb{S}$, a.s.. Furthermore, for any $p > 0$, there exist constants $K_{i}(p) > 0$ such that

Proof. It is easy to see the function defined in right side of the system (1) obeys the local Lipschitz condition. Then, the system has a unique local solution $X(t)$ on $[0, \tau_{e})$, where $\tau_{e}$ stands for the explosion time. We may prove $\tau_{e} = \infty$ a.s.. Fix a $k_{0} > 0$ sufficiently large for $x_{1}(\phi), x_{2}(\phi), x_{3}(\phi)\in (1/k_{0}, k_{0})$. For each integer $k > k_{0}$, define stopping times as follows

It is clear that $\tau_{k}$ is increasing with $k$. Setting $\tau_{\infty} = \lim\limits_{t \to \infty}\tau_{k}$, we have $\tau_{\infty}\le \tau_{e}$ a.s.. Thus, we only need to prove $\tau_{\infty} = \infty$, a.s.. For if this statement is false, there exists a $T > 0$ and an $\epsilon\in (0, 1)$ such that $P(\tau_{\infty}\le T) > \epsilon$. We can find an integer $k_{1} > k_{0}$ such that $P(\tau_{k}\le T) > \epsilon$ for any $k > k_{1}$. Then take three positive constants $\chi, \kappa$ and $n > 0$ for

We choose $n$ sufficiently large and two positive constants $\kappa > \frac{2c_{22}}{c_{32}}$ and $\chi > \frac{nc_{21}+n\kappa c_{31}}{2c_{11}}$. Define $v_{i}(x_{i}) = x_{i}-1-\ln x_{i}, i = 1, 2, 3,$

In particular, $v_{i}(x_{i}), i = 1, 2, 3$ and $v_{4}(x_{1}, x_{2}, x_{3})$ are independent of $\gamma$. By the generalized Itô's formula, we can obtain

where

Thus

Define $V(x_{1}, x_{2}, x_{3}) = \chi v_{1}(x_{1})+v_{2}(x_{2})+\kappa v_{3}(x_{3})+v_{4}(x_{1}, x_{2}, x_{3})$

where

Thus there is a $L > 0$

Integral on both sides of Eq (3) from $0$ to $\tau_{k} \wedge T$ and then take expectation, for $k\rightarrow \infty$, we can get the following contradiction

Thus, $\tau_{e} = \infty$, a.s.. We next to prove $\limsup\limits_{t \to \infty} E[x_{i}^{p}(t)]\leq K_{i}(p), i = 1, 2, 3.$ Define a function $U_{1}(x_{1}) = e^tx_{1}^{p}.$ According to the generalised Itô's formula

where

For $p\ge 1$

where $K^*_{1}(p) = \sup\limits_{x_{1} > 0}\left\{\left[1+\frac{p(p-1)\hat \sigma_{1}^{2}}{2}+p\left(\hat a_{1}-h_{1}\right)\right]x_{1}^{p}-pc_{11}x_{1}^{p+1}\right\}$ is a constant. For $0 < p < 1$

where $K^{**}_{1}(p) = \sup\limits_{x_{1} > 0}\left\{\left(1+p\hat a_{1}\right)x_{1}^{p}-pc_{11}x_{1}^{p+1}\right\}$. Integral on both sides of the Eq (4) from $0$ to $t$ and then take the expectation, We can then show that $E[e^{t}x_{1}^{p}(t)]\leq x_{1}^{p}(0)+K_{1}(p)(e^{t}-1),$ where $K_{1}(p) = \max\{K^*_1(p), K^{**}_1(p)\}.$ Thus, $\limsup\limits_{t\to \infty}E[x_{1}^{p}(t)]\leq K_{1}(p).$ Continuing this approach we define

Letting an appropriate $n > 0$ such that $c_{22}-\frac{pc_{21}}{1+p}n^{-\frac{1+p}{p}} > 0, p(c_{33}-\frac{pc_{31}}{1+p}n^{-\frac{1+p}{p}}-\frac{pc_{32}}{1+p}n^{-\frac{1+p}{p}}) > 0.$ The constants $c_{1}^{*}, {c_{2}}^{*}$ and ${c_{3}}^{*}$ satisfying

By the generalized Itô's formula, we obtain

where

For $p > 1$

where $K^*_{2}(p) = \sup\limits_{x_{1}, x_2 > 0}\Bigg\{c_{1}^{*}\Bigg[1+\frac{p(p-1){\hat \sigma_{1}}^{2}}{2}+p\hat{a_{1}}\Bigg]x_{1}^{p}+\Bigg[1+\frac{p(p-1){\hat \sigma_{2}}^{2}}{2}\Bigg]x_{2}^{p}-\frac{p^{2}e^{\tau_{21}}}{1+p}n^{1+p}c_{21}x_{1}^{1+p}-p(c_{22}-\frac{pc_{21}}{1+p}n^{-\frac{1+p}{p}})x_{2}^{1+p}\Bigg\}$, and for $0 < p < 1$

where $K^{**}_2(p) = \sup\limits_{x_{1}, x_2 > 0}\Bigg\{c_{1}^{*}\Big(1+p\hat{a_{1}}\Big)x_{1}^{p}+x_{2}^{p}-\frac{p^{2}e^{\tau_{21}}}{1+p}n^{1+p}c_{21}x_{1}^{1+p}-p(c_{22}-\frac{pc_{21}}{1+p}n^{-\frac{1+p}{p}})x_{2}^{1+p}\Bigg\}$. Take the same method as above, we obtain there's a positive constant $K_{2}(p) = \max\{K^*_2(p), K^{**}_2(p)\}$ such that $\limsup\limits_{t\to \infty} E[x_{2}^{p}(t)]\leq K_{2}(p).$ Next, we define

Then

where

For $p > 1$

For $0 < p < 1$

Similarly, there's a positive constant $K_{3}(p)$ satisfying $\mathcal{L}U_{3}(t)\leq K_{3}(p)e^{t}$ and $\limsup\limits_{t\to \infty} E[{x_{3}}^{p}(t)]\leq K_{3}(p).$ This completes the proof.

In the following, we give some basic lemmas.

Lemma 2. ^[41] Let $Z(t)\in \mathbb{R}_{+}$ and $g(t)$ be two stochastic processes satisfying $\lim\limits_{t \to \infty}\frac{g(t)}{t} = 0$, a.s.. We have the following conclusion.

(ⅰ) If there are three constants $T > 0, l_{2} > 0$ and $l_{1}$ such that for all $t\ge T$

then

(ⅱ) If there exist three positive constants $T, l_{1}$and $l_{2}$ such that

for all $t\ge T$, then

Consider the following auxiliary system:

with initial date $\varphi \in \mathbb{C}\left([-\tau, 0], \mathbb {R}_{+}^{3}\right)$ with $\gamma(0) = i \in \mathbb{S}$. Obviously, the system (5) has a unique global positive solution ^[23].

Lemma 3. For the system (5), we have the following conclusion.

(a) If ${\overline u }_1 < 0$, then $\lim\limits_{t \to \infty}Y_{j}(t) = 0, a.s., j = 1, 2, 3;$

(b) If ${\overline u }_1 = 0$, then $\langle Y_{1}(t) \rangle = 0, \lim\limits_{t \to \infty}Y_{i}(t) = 0, a.s., i = 2, 3;$

(c) If ${\overline u }_1 > 0, \Delta_2-\tilde{\Delta}_2 < 0$, and $\Delta_3-\tilde{\Delta}_3 < 0$, then

(d) If ${\overline u }_1 > 0, \Delta_2-\tilde{\Delta}_2 < 0$ and $\Delta_3-\tilde{\Delta}_3\ge 0$, then

(e) If ${\overline u }_1 > 0, \Delta_2-\tilde{\Delta}_2\ge 0$ and $c_{22}(\Delta_3-\tilde{\Delta}_3)+c_{32}(\Delta_2-\tilde{\Delta}_2) < 0$, then

(f) If ${\overline u }_1 > 0, \Delta_2-\tilde{\Delta}_2\ge 0$ and $c_{22}(\Delta_3-\tilde{\Delta}_3)+c_{32}(\Delta_2-\tilde{\Delta}_2)\ge 0$, then

Proof. Firstly, let us prove the conclusion $(a)$. By the generalized Itô's formula, we have

Because of the ergodicity of $\gamma(t)$, one gets

Obviously, for any $\epsilon \in (0, \bar u_{1})$, there exist a set $\Omega_{\epsilon}\subset \Omega$ satisfying $P(\Omega_{\epsilon})\ge 1-\epsilon$ and a positive constant $T = T(\epsilon)$ such that for $t\ge T$

Substituting Eq (10) into Eq (6), then

Meanwhile, the quadratic variation of the stochastic integral $\int^{t}_{0}\sigma_{j} \left(\gamma\left(s\right)\right)\mathrm{d}\beta_{j}(s)$ is $\int^{t}_{0}\sigma_{j}^{2} \left(\gamma\left(s\right)\right)\mathrm{d}s \leq {\hat\sigma_{j}}^{2}t$. The strong law of large number theorem shows

When ${\overline u }_1 < 0$, for Eq (12) and Lemma 2, we can show that $\lim\limits_{t\to \infty}Y_{1}(t) = 0, \lim\limits_{t \to \infty} t^{-1}\int_{0}^{t}Y_{1}(s)\mathrm{d}s = 0$. From Eqs (7) and (8), for $t \to \infty$, then

For any $\epsilon \in (0, \bar u_{2})$ satisfying $(-\bar u_{2}+\epsilon) < 0$, by Eqs (9) and (13), we have

It follows from Lemma 2 that $\lim\limits_{t \to \infty}Y_{2}(t) = 0$ and $\langle Y_{2}(t) \rangle = 0$, for any $t\ge T$. From Eq (14),

It's easy to see $\lim\limits_{t \to \infty}Y_{3}(t) = 0, a.s.$. This completes the proof of (a).

Now we are in the position to prove (b). By Eq (11) and Lemma 2, we have

According to the arbitraryiness of $\epsilon$, we get $\langle Y_{1}(t) \rangle = 0$ when ${\overline u }_1 = 0$. Similarly, we obtain the Eqs (13) and (14) from Eqs (7) and (8). This implies $\lim\limits_{t \to \infty}Y_{j}(t) = 0, j = 2, 3$. The proof of (b) is complete.

Next, we shall prove (c). It can be shown from Eq (12) that $\langle Y_{1}(t) \rangle = \frac{{\overline u }_1}{c_{11}}\; \ a.s.$ when ${\overline u }_1 > 0$. Fix a positive constant $T$, for any $t \ge T$, we may shift Eq (8) to obtain

where $\phi_{2}(t) = c_{21}\int^{0}_{-\tau_{21}}\int^{0}_{\theta}Y_{1}(s)\mathrm{d}s\mathrm{d}\mu_{21}(\theta)-c_{21}\int^{0}_{-\tau_{21}}\int^{t}_{t+\theta}Y_{1}(s)\mathrm{d}s\mathrm{d}\mu_{21}(\theta)$ and then $\lim\limits_{t \to \infty}t^{-1}\phi_2(t) = 0$. When $\Delta_2-\tilde{\Delta}_2 < 0$, we further have $\lim\limits_{t \to \infty}Y_{2}(t) = 0, a.s.$ by Lemma 2. For any $t\ge T$, Eq (8) follows

where $\phi_{3}(t) = c_{31}\int^{0}_{-\tau_{31}}\int^{0}_{\theta}Y_{1}(s)\mathrm{d}s\mathrm{d}\mu_{31}(\theta)-c_{31}\int^{0}_{-\tau_{31}}\int^{t}_{t+\theta}Y_{1}(s)\mathrm{d}s\mathrm{d}\mu_{31}(\theta).$ Clearly, $\lim\limits_{t \to \infty}\frac{1}{t}\phi_{3}(t) = 0$. From Lemma 2, when $\Delta_3-\tilde{\Delta}_3 < 0, \lim\limits_{t \to \infty}Y_{3}(t) = 0, a.s.$. Thus the required conclusion (c) follows.

Next, we shall give the proof of (d). Conclusion (c) shows $\langle Y_{1}(t) \rangle = \frac{{\overline u }_1}{c_{11}}$ and $\lim\limits_{t \to \infty}Y_{2}(t) = 0, a.s.$ for ${\overline u } > 0$ and $\Delta_2-\tilde{\Delta}_2 < 0$. Assume $\Delta_3-\tilde{\Delta}_3 \ge 0$, there is a conclusion that $\langle Y_{3}(t) \rangle = \frac{ \Delta_3-\tilde{\Delta}_3}{c_{11}c_{33}}$ from Eq (16). The proof of (d) is completed.

In the following, we shall prove (e). From Eqs (11) and (15), assume $\Delta_2-\tilde{\Delta}_2 \ge 0$, we get $\langle Y_{1}(t) \rangle = \frac{{\overline u }_1}{c_{11}}$ and $\langle Y_{2}(t) \rangle = \frac{\Delta_2-\tilde{\Delta}_2}{c_{11}c_{22}}$. Then, for any $t \ge T$, we have

where

Then, $\lim\limits_{t \to \infty}t^{-1}\phi^{*}_3(t) = 0$. When $c_{22}(\Delta_3-\tilde{\Delta}_3)+c_{32}(\Delta_2-\tilde{\Delta}_2) < 0$, we get $\lim\limits_{t \to \infty}Y_{3}(t) = 0$. This completes the proof of (e).

Lastly, we give the proof of (f). Let ${\overline u } > 0$ and $\Delta_2-\tilde{\Delta}_2 > 0$, we can then show that $\langle Y_{1}(t) \rangle = \frac{{\overline u }_1}{c_{11}}$ and $\langle Y_{2}(t) \rangle = \frac{\Delta_2-\tilde{\Delta}_2}{c_{11}c_{22}}\; \ a.s.$ by Eqs (11) and (15). From Eq (17), we directly obtain $\langle Y_{3}(t) \rangle = \frac{c_{22}\Gamma_{3}+c_{32}\Gamma_{2}}{c_{11}c_{22}c_{33}}$ provided that $c_{22}(\Delta_3-\tilde{\Delta}_3)+c_{32}(\Delta_2-\tilde{\Delta}_2)\ge 0$.

This completes the proof of all cases.

Lemma 4. Let $X(t) = \left(x(t), \gamma(t)\right)\in \mathbb{R}^{3}_{+} \times \mathbb{S}$ be a global positive solution of the system (1). Then the solution has the following properties

Proof. It is easy to see

By comparison theorem, we obtain $x_{j}(t) \leq Y_{j}(t), \; a.s., j = 1, 2, 3$. It follows from Lemma 3 that $\lim\limits_{t \to \infty}Y_{j}(t) = 0$ or $\langle Y_{j}(t)\rangle = a, j = 1, 2, 3$, where the constant $a\geq 0$. From Eqs (6)–(8), then

This completes the proof.

Lemma 5. Let condition ${\bf (H_{2})}$ hold. For the defined parameters $R, \omega_i, \Theta_i, i = 1, 2, 3$, we have the following conclusions.

(a). If $R > 0$, then $\omega_1 > \omega_2 > \omega_3$.

(b). If $R < 0$, then $\omega_1 > \Theta_2 > \Theta_3$.

(c). If $R = 0$, then $\omega_1 > \Theta_2 = \Theta_3$.

Proof. From the definition of $R, \omega_i, \Theta_i, i = 1, 2, 3$, we have

When $R > 0$, we can see $\omega_1 > \omega_2 > \Theta_2 > \omega_3$. When $R < 0$, then $\Theta_3 < \omega_2 < \Theta_2 < \omega_1$. When $R = 0$, then $\omega_3 = \Theta_3 = \omega_2 = \Theta_2 < \omega_1$. Thus the proof is complete.

3.   Global analysis for dynamical behaviors

In this section, we shall investigate dynamical behaviors involving persistence, extinction, stability and attractiveness of the system (1). In terms of parameters $R$, $\omega_i$ and $\Theta_i$, we first give a global analysis of dynamical behaviors according to Lemma 5.

Theorem 1.\quad Let conditions ${\bf (H_{1})}$ and ${\bf (H_{2})}$ be satisfied and $X(t) = \left(x(t), \gamma(t)\right)\in \mathbb{R}^{3}_{+} \times \mathbb{S}$ be a global positive solution of the system (1). We have the following conclusions.

(ⅰ). When $R > 0$

(a). If $\omega_1 < 1$, then $\lim\limits_{t \to \infty}x_{i}(t) = 0, \; a.s., \; i = 1, 2, 3;$

(b). If $\omega_1 = 1$, then $\langle x_{1}(t)\rangle = 0, \lim\limits_{t \to \infty}x_{j}(t) = 0, a.s., j = 2, 3;$

(c). If $\omega_1 > 1 > \omega_2$, then

(d). If $\omega_2 = 1$, then

(e). If $\omega_2 > 1 > \omega_3$, then

(f). If $\omega_3 = 1$, then

(g). If $\omega_3 > 1$, then

(ⅱ). When $R < 0$

(h). If $\omega_1 < 1$, then $\lim\limits_{t \to \infty}x_{i}(t) = 0, \; a.s., \; i = 1, 2, 3;$

(i). If $\omega_1 = 1$, then $\langle x_{1}(t)\rangle = 0, \lim\limits_{t \to \infty}x_{j}(t) = 0, a.s., j = 2, 3;$

(j). If $\omega_1 > 1 > \Theta_2$, then

(k). If $\Theta_2 = 1$, then

(l). If $\Theta_2 > 1 > \Theta_3$, then

(m). If $\Theta_3 = 1$, then

(n). If $\Theta_3 > 1$, then

(ⅲ). When $R = 0$

(o). If $\omega_1 < 1$, then $\lim\limits_{t \to \infty}x_{i}(t) = 0, \; a.s., \; i = 1, 2, 3;$\\

(p). If $\omega_1 = 1$, then $\langle x_{1}(t)\rangle = 0, \lim\limits_{t \to \infty}x_{j}(t) = 0, \; a.s., j = 2, 3;$\\

(q). If $\Theta_2 < 1 < \omega_1$, then

(r). If $\Theta_2 = 1$, then

(s). If $\Theta_2 > 1$, then

Proof. By the stochastic comparison theorem and the condition (a) in Lemma 3. when $\omega_1 < 1$, we have the conclusion

Next, we prove the conclusion (b) holds. Thanks to the condition (b) in Lemma 3, we can immediately obtain that $\langle x_{1}(t)\rangle \leq \langle Y_{1}(t) \rangle = 0$ and $\lim\limits_{t \to \infty}x_{i}(t) \leq \lim\limits_{t \to \infty}Y_{i}(t) = 0, a.s., i = 2, 3.$

Then, we shall prove the condition (c). When $\omega_1 > 1 > \omega_2$, that is, $\bar{u}_1 > 0, \Delta_2-\tilde{\Delta}_2 < 0$ and $\Delta_3-\tilde{\Delta}_3 < 0$. we further get $\lim\limits_{t \to \infty}x_{j}(t) \leq \lim\limits_{t \to \infty}Y_{j}(t) = 0, a.s., j = 2, 3$ by condition (c) in Lemma 3. Fix a constant $T > 0$ and any $t \ge T$, Applying generalized Itô's formula to the system (1) leads to

where

When $\lim\limits_{t \to \infty}x_{j}(t), j = 2, 3$, for Eq (18), one can observe that

From Lemma 2, we have $\langle x_{1}\rangle = \frac{\bar{u}_1}{c_{11}}$. This completes the proof of (c).

Then, we shall prove the conclusion (d). Similar to the above method, when $\omega_2 = 1$, that is, $\bar{u}_1 > 0, \Delta_2-\tilde{\Delta}_2 = 0$ and $\Delta_3-\tilde{\Delta}_3 < 0$. The condition (e) in Lemma 3 shows $\lim\limits_{t \to \infty}Y_{3}(t) = 0$ and $\langle Y_2(t) \rangle = 0$. Thus we can check that $\langle x_2(t) \rangle \leq \langle Y_2(t) \rangle = 0$ and $\lim\limits_{t \to \infty}x_{3}(t) \leq \lim\limits_{t \to \infty}Y_{3}(t) = 0$. Eq (18) changes into Eq (21) then we get $\langle x_{1}(t)\rangle = \frac{\bar{u}_1}{c_{11}}$.

Next, we prove the conclusion (e). It is not difficult to show that $\lim\limits_{t \to \infty}t^{-1}\Phi_{i}(t) = 0, i = 1, 2, 3$. Take two negative constants $m$ and $n$ satisfying

Then

Compute $m\times$ Eq (18) + Eq (19) + $n$ $\times$ Eq (20)

From Lemma 4, there exists an $\epsilon > 0$ and a $T = T(\epsilon) > 0$. For $t \ge T$, one has that

Because of $\frac{C_2}{\tilde{C}_2} > \omega_2 > 1$, it's easy to see $C_2-\tilde{C}_2 > 0$. we further have $\sup \langle x_{2}(t) \rangle \leq \frac{C_2-\tilde{C}_2}{|C|}$. Then take two negative constants $m^*$ and $n^*$ satisfying

That is

Compute Eq (18) + $m^*\times$ Eq (19) + $n^*\times$ Eq (20)

Thus

From $({\bf H_{2}})$, we gain $C_1-\tilde{C}_1 > 0$ and then $\sup \langle x_{1}(t) \rangle \leq \frac{C_1-\tilde{C}_1}{|C|}.$ For Eq (20)

Assume $\omega_3 = \frac{C_3}{\tilde{C}_3} < 1$, we obtain $\lim\limits_{t \to \infty}x_{3}(t) = 0$. Equations (19) and (20) follow that:

where

Analogously, $\lim\limits_{t \to \infty}t^{-1}\Phi_i^*(t) = 0, a.s., i = 1, 2.$ Compute $c_{21} \times$ Eq (30) + $c_{11} \times$ Eq (31)

For any $t \ge T$

It means $\sup \langle x_{2}(t) \rangle \leq \frac{\Delta_2-\tilde{\Delta}_2}{C_{33}}$ and then

we can see if $\omega_1 > 1$,

Then $\inf \langle x_{1}(t) \rangle \ge \frac{\Delta_1-\tilde{\Delta}_1}{C_{33}}$ and

We can now easily establish $\langle x_{2}(t) \rangle = \frac{\Delta_2-\tilde{\Delta}_2}{C_{33}}$ and then

Finally, we have $\langle x_{1}(t) \rangle = \frac{\Delta_1-\tilde{\Delta}_1}{C_{33}}, \langle x_{2}(t) \rangle = \frac{\Delta_2-\tilde{\Delta}_2}{C_{33}}$ and $\lim\limits_{t \to \infty}x_{3}(t) = 0.$ This completes the proof of (e).

Next, we prove the case (f). For Eq (24), If $\omega_3 = 1$, this is, $C_3-\tilde{C}_3 = 0$, then we get $\sup\limits_{t \to \infty}\langle x_3(t) \rangle \leq 0.$ Because of $\langle x_3(t) \rangle \ge 0$, we finally derive $\langle x_3(t) \rangle = 0$. Equations (18) and (19) become Eqs (25) and (26) respectively. Finally we derive $\langle x_{1}(t) \rangle = \frac{\Delta_1-\tilde{\Delta}_1}{C_{33}}, \langle x_{2}(t) \rangle = \frac{\Delta_2-\tilde{\Delta}_2}{C_{33}}$ and $\langle x_3(t) \rangle = 0$.

Lastly, we shall prove the condition (g) holds. If $\omega_3 = \frac{C_3}{\tilde{C}_3} > 1$, for Eq (29), we have $\sup \langle x_{3}(t) \rangle \leq \frac{C_3-\tilde{C}_3}{|C|}.$ and

Thus, $\frac{C_1-\tilde{C}_1}{|C|} \leq \inf \langle x_{1}(t)\rangle \leq \sup \langle x_{1}(t)\rangle \leq \frac{C_1-\tilde{C}_1}{|C|}$, $\langle x_{1}(t)\rangle = \frac{C_1-\tilde{C}_1}{|C|}$. By the same way, for Eq (20)

We have $\frac{C_2-\tilde{C}_2}{|C|} \leq \inf \langle x_{2}(t)\rangle \leq \sup \langle x_{2}(t)\rangle \leq \frac{C_2-\tilde{C}_2}{|C|}$ and then $\langle x_{2}(t)\rangle = \frac{C_2-\tilde{C}_2}{|C|}$. Similarly, $\langle x_{3}(t)\rangle = \frac{C_3-\tilde{C}_3}{|C|}$.

The proof of (ⅱ) and (ⅲ) are similar to one of (ⅰ). Here we omit the remaining part and all the proof is complete. In the following, we further give the attractiveness of all positive solutions by $M$-matrix.

Theorem 2. If the following matrix is a non-singular $M$-matrix

then

where $X^{\eta, k}(t)$ and $X^{\eta^*, k}(t)$ are two positive solutions of the system (1) with initial conditions $\eta \in \mathbb{C}\left([-\tau, 0], \mathbb {R}_{+}^{3}\right), \gamma(0) = k$ and $\eta^* \in \mathbb{C}\left([-\tau, 0], \mathbb {R}_{+}^{3}\right), \gamma(0) = k$.

Proof. Notice that

Since $C_M$ is an non-singular $M$-matrix, there exists positive vector $\zeta = (\zeta_1, \zeta_2, \zeta_3)^T, \zeta_{i} > 0, i = 1, 2, 3$ such that $C_M \zeta > 0$. We just need to proof $\lim\limits_{t \to \infty}E|(x_{j}(t; \eta)-x_{j}(t; {\eta}^{*})| = 0, j = 1, 2, 3,$ define

Calculating the upper right derivative of $v_{j}(x_{j})$ along the solution of Eq (33), it follows that

Let

Then

Integrate both sides of the Eq (29) from $0$ to $t$ and then take expectation

From $C_M\zeta > 0$, we see $\int^{t}_{0}E|x_{j}(s, \eta)-x_{j}(s, {\eta}^{*})| \mathrm{d}s < \infty, j = 1, 2, 3.$ Next, we define function

For arbitrary $t_{1}, t_{2} \in [-\tau, +\infty)$,

Integrating both sides of the system (1) from $t_{1}$ to $t_{2}$, we see

For any $t_{2} > t_{1}, p > 0$

Moreover,

where

Similarly,

In view of the Theorem 2.11 in ^[40], for any $t_{2} > t_{1}$ and $p\ge 2$, we have

From Lemma 1 and Eqs (30)–(34), for $p \ge 2$, we get

where, for $|t_2-t_1| < \epsilon$, the constant

By the same method, we get ${\left[E|x_{1}(t_{2}, \eta^*)-x_{1}(t_{1}, \eta^*)|\right]}^{p} \leq L_{1}|t_{2}-t_{1}|^{\frac{p}{2}}.$ Thus $|G_{1}(t_{2})-G_{1}(t_{1})| \leq 2{L_{1}}^\frac{1}{p}\sqrt{|t_{2}-t_{1}|}.$ Similarly, there exist two constants $L_{2} > 0$ and $L_{3} > 0$ such that

So the $G_{j}(t), j = 1, 2, 3,$ is uniformly continuous function. According to Barb\v{a}lat Lemma ^[42], we have the conclusion. The proof is complete.

4.   The optimal harvesting policy

In this section, we shall provide the sufficient and necessary conditions ensuring the existence of optimal solution of the system (1) with harvesting.

Assume $y(t)$ is the stochastic process of $(X(t), \gamma(t))$ and $p(t, \phi, i, \mathrm{d}y \times \{j\})$ is the transition probability of the process $\gamma(t)$. $P(t, \phi, i, A \times D)$ stands for the probability of event $\{y(t) \in A \times D\}$ with the initial condition $\phi \in \mathbb{C}\left([-\tau, 0], \mathbb {R}_{+}^{3}\right), \gamma(0) = i$ and

Now, let $\mathcal{P}$ stands for all probability measures on $\mathbb{C}\left([-\tau, 0], \mathbb {R}_{+}^{3}\right) \times \mathbb{S}$. For any measures $P_{1}, P_{2} \in \mathcal{P}$ define the metric

where

Theorem 3.\quad Let all positive solutions of (1) be globally attractive. Then the solution $(X(t), \gamma(t))$ has a unique ergodic invariant distribution $\bar{v}(\cdot \times \cdot)$ defined on $\mathbb{C}\left([-\tau, 0], \mathbb {R}_{+}^{3}\right) \times \mathbb{S}$.

Proof. We observe from Lemma 1 that for any initial data $\phi \in \mathbb{C}\left([-\tau, 0], \mathbb {R}_{+}^{3}\right)$ with $\gamma(0) = i \in \mathbb{S}$, the family of transition probability $\{p(t, \phi, i, \mathrm{d}y \times \{j\}): t\ge 0\}$ is tight by the Chebyshev inequality. That is to say, for a compact subset $K = K(\epsilon, \phi, i)$ and any $\epsilon > 0$

Next, based on the Theorem 2, for any $\epsilon > 0, t\ge T$ and any compact subset $K \in \mathbb{C}\left([-\tau, 0], \mathbb {R}_{+}^{3}\right)$, we have

uniformly in $\phi, \tilde{\phi} \in K$ and $i, j \in \mathbb{S}$. The proof is similar as Lemma 5.6 in ^[40], here we omit it. And then for an arbitrary initial data $\eta \in \mathbb{C}\left([-\tau, 0], \mathbb {R}_{+}^{3}\right)$ with $\gamma(0) = j \in \mathbb{S}$, the $\{p(t, \eta, j, \cdot \times \cdot):t\ge 0\}$ is cauchy in the metric space $\mathcal{P}(\mathbb{C}\left([-\tau, 0], \mathbb {R}_{+}^{3}\right) \times \mathbb{S})$ (see Lemma 5.7 in ^[40]). Hence, there exists a unique probability measure $\bar{v}(\cdot)$ such that

The last, for any initial value $\phi \in \mathbb{C}\left([-\tau, 0], \mathbb {R}_{+}^{3}\right)$ with $\gamma(0) = i \in \mathbb{S}$, we obtain that

It means that the $p(t, \phi, i, \cdot \times \cdot)$ converges weakly to $\bar{v}(\cdot)$. ^[43] (Corollary 3.4.3 and Theorem 3.2.6) show $\bar{v}(\cdot)$ is strong mixing and ergodic.

This completes the proof.

When $C_M$ is a non-singular matrix, it follows from Theorems 2 and 3 that any positive solution $(X(t), \gamma(t))$ of the system (1) has a unique ergodic invariant distribution $\bar{v}(\cdot \times \cdot)$ defined on $\mathbb{C}\left([-\tau, 0], \mathbb {R}_{+}^{3}\right)\times \mathbb{S}$. This guarantees that $\lim\limits_{t \to \infty}\sum\limits^{3}_{j = 1}E(h_{j}x_{j}(t))$ exists for the positive solution $x(t) = (x_1(t), x_2(t), x_3(t))^T$. That is, the harvesting effort $H^{*} = (h_{1}^{*}, h_{2}^{*}, h_{3}^{*})^T$ is the optimal solution of Eq (2) if

$(i) \Phi(H^*)$ is the maximum value of $\Phi$;

$(ii)$ All species in the system (1) are stable in mean (see also Definition 1).

The following theorem gives sufficient and necessary conditions ensuring the existence of optimal solution $H^* = (h_{1}^{*}, h_{2}^{*}, h_{3}^{*})^T$ for Eq (2).

Theorem 4. Let $C_M$ be a non-singular matrix. Then the system (1) has the optimal harvesting policy if and only if (a), (b) and (c) all hold.

Furthermore, we obtain the optimal harvesting effort $H^{*} = \left[\left(C(C^{-1})^{T}+I\right)\tilde{R}\right]^{-1}d$ and the maximum of sustainable yield $\Phi(H^{*}) = (H^*)^{T}\tilde{R}C^{-1}(d-\tilde{R}H^*)-W$, where $\tilde{R} = diag\{r_1, r_2, r_3\}$.

Proof. Denote $\mu = \{H = (h_{1}, h_{2}, h_{3})^{T} \in \mathbb R^3\mid (a), (b)\; \ \text{all hold}\}$. If we take $h_1 = h_2 = h_3 = 0, r_1 = r_2 = r_3 = 1$ and choose the appropriate parameters, we can see $\mu$ is not empty. Denote $\rho$ means the stationary probability density of model (1) and $x(t) = (x_1(t), x_2(t), x_3(t))^T$. For notation convenience, we use $x$ instead of $x(t)$, so we get

Theorem 3 shows the system (1) has a unique stable distribution $\bar{v}$. Due to the one-to-one corresponded between $\rho$ and $\bar{v}$, we have the following relations (see Theorem 3.3.1 in ^[43])

For arbitrary $(h_{1}, h_{2}, h_{3})^{T} \in \mu$, we have

Thus

Next

Assume $H^* = (h^*_{1}, h^*_{2}, h^*_{3})^{T}$ is the solution of

Thus one get

The Hessian matrix

Due to $C^{-1}+(C^{-1})^{T}$ is positive definite, then $H^*$ is a unique extreme point of $\Phi(H)$. And then the optimal harvesting effort $H^* = \left[\left(C(C^{-1})^{T}+I\right)\tilde{R}\right]^{-1}d$ and the maximum value of expectation of sustainable yield $\Phi(H^{*}) = (H^*)^{T}\tilde{R}C^{-1}(d-\tilde{R}H^*)-W$.

If conditions (a) and (b) do not hold, the species of the system (1) will not stable in mean or to be extinct. Then the optimal harvesting policy does not exist. What's more, if (c) not hold. That is to say $C^{-1}+(C^{-1})^{T}$ is not positive definite. Set $Q = C^{-1}+(C^{-1})^{T}$, then we can see

So $Q$ must not be negative definite, and $Q$ is indefinite. This implies that $\Phi(H)$ has no extreme point. Thus, the optimal harvesting policy does not exist. The proof is complete.

5.   Examples and numerical simulations

In this section, we shall give some numerical examples and their simulations to illustrate the effectiveness of the obtained results for the hybrid stochastic phytoplankton-zooplankton-fish model with distributed delays (1). For the practical biological significance, some parameter values of the system (1) are given such as the growth rate of phytoplankton, the death rate of zooplankton and fish ^[44,45,46]. According to these publications, we take the mean growth rate of phytoplankton ranged from $0.26d^{-1}$ to $1.04d^{-1}$, the death rate of diaphanosoma brachyurum (a kind of zooplankton) ranged from 0.01$d^{-1}$ to 0.92$d^{-1}$, and the death rate of fish ranged from 0.15$d^{-1}$ to 0.95$d^{-1}$. Due to technical reasons, the rate of intra-specific competition, capture and food conversation are not easy to be obtained. In the following examples, we will refer to the value range of system parameters.

Example 1. Consider the stochastic delayed phytoplankton-zooplankton-fish model (1) with $\mathbb{S} = \{1, 2\}, \xi = (0.3, 0.7)$, and parameters $a_{1}(1) = 0.56, a_{2}(1) = 0.27, a_{3}(1) = 0.45$, $a_{1}(2) = 0.80, a_{2}(2) = 0.33, a_{3}(2) = 0.64$, $h_{1} = h_{2} = h_{3} = 0$, $c_{11} = 0.35, c_{12} = 0.2, c_{13} = 0.1, c_{21} = 0.15, c_{22} = 0.65, c_{23} = 0.3, c_{31} = 0.28, c_{32} = 0.1, c_{33} = 0.55$, $\sigma_{1}^2(1) = \sigma_{1}^2(2) = 0.02, \sigma_{2}^2(1) = \sigma_{2}^2(2) = 0.01, \sigma_{3}^2(1) = \sigma_{3}^2(2) = 0.012$.

We compute that $C_{23} = 0.021 > 0, C_{31} = -0.005 < 0, C_{12} = 0.0015 > 0, |C| = 0.155025 > 0$, $\bar{a}_{1} = 0.73, \bar{a}_{2} = 0.31, \bar{a}_{3} = 0.58, C_1 = 0.323, C_2 = 0.00015, C_3 = 0.00097, \tilde{C}_1 = 0.0032, \tilde{C}_2 = 0.000367, \tilde{C}_3 = 0.0034, R = 0.000001 > 0, \omega_3 = 0.2845 < 1, \omega_2 = 0.3077 < 1, \omega_1 = 73 > 1$. It follows from Case (c) of Theorem 1 that

Figure 2 is shown to illustrate the species $x_2, x_3$ are extinct and $x_1$ is stable in mean. So, there is no optimal solution for the optimization problem (2).

The above example shows the effect of time delay and environmental noise on the survival and extinction of species without harvesting. According to Theorem 1, the time delay has no effect on the survival and extinction of the three species, but the harvesting and environmental noise have great influence on it. This is also verified by our numerical experiments. The results of example 1 showed that when Zooplankton and Fish were extinct, Phytoplankton could still stable in mean under certain environmental noise. But it is meaningless to discuss the optimal harvesting policy in the case of species extinction.

In order to demonstrate some results about optimal harvesting policy for the system (1), we give the following example.

Example 2. uad Consider the stochastic delayed phytoplankton-zooplankton-fish model (1) with $\mathbb{S} = \{1, 2\}, \xi = (0.9, 0.1)$, $a_{1}(1) = 0.9789, a_{2}(1) = 0.0211, a_{3}(1) = 0.1044$, $a_{1}(2) = 0.59, a_{2}(2) = 0.01, a_{3}(2) = 0.06$, we compute $\bar{a}_{1} = 0.94, \bar{a}_{2} = 0.02, \bar{a}_{3} = 0.1$.

When there is no harvesting, i.e. $h_i = 0, i = 1, 2, 3$, we compute that $C_1 = 0.3672, C_2 = 0.0062, C_3 = 0.1599, \tilde{C}_1 = 0.0032, \tilde{C}_2 = 0.000367, \tilde{C}_3 = 0.0034, R = -0.000243 < 0, \Theta_3 = \frac{C_2}{\tilde{C}_2} = 16.8163 > 1$. It follows from Case (n) of Theorem 1 that

Figure 3 is shown to illustrate the species $x_i, i = 1, 2, 3$ are stable in mean for the system (1) with no harvesting.

Furthermore, we select the appropriate values of the unit profit parameters $r_1 = 0.6, r_2 = 2.7, r_3 = 4$, and fixed cost $W = 0.2$. Then, all the conditions in Theorem 4 are checked as follows: $d = (0.93, 0.025, 0.106)^{T}, H^* = (h^*_{1}, h^*_{2}, h^*_{3})^{T} = \left[\left(C(C^{-1})^{T}+I\right)\tilde{R}\right]^{-1}d = (0.6915, 0.0713, 0.0764)^{T}, \tilde{C}_1|_{H = H^*} = 0.2622, \tilde{C}_2|_{H = H^*} = 0.0059, \tilde{C}_3|_{H = H^*} = 0.1578, R|_{H = H^*} = 0.000225 > 0, \omega_3 = \frac{C_3}{\tilde{C}_3} = 1.0129 > 1$ and $C_M^{-1}\geq 0$. Thus, the optimal harvesting effort is $H^{*} = (0.6915, 0.0713, 0.0764)^T$ and the maximum of total economic income is 0.4792. Under the optimal harvesting condition, by Case (g) in Theorem 1, we see that

Figure 4 is shown to illustrate the species $x_i, i = 1, 2, 3$ are stable in mean for the system (1) with the optimal harvesting $H^*$.

In example 2, we see that all species are stable in mean without harvesting. According to Theorem 4, we get the maximum harvesting effort $H^*$, and verify that all species are still stable in mean. This means that the system (1) has an optimal harvesting strategy with the maximum of total economic income. However, excessive harvesting or strong enough environmental noise may cause the extinction of species, and the optimal harvesting strategy does not exist.

Acknowledgments

This work is partially supported by the National Natural Science Foundation of China (NSFC) under Grant No.11971081, No. 11701060, the Fundamental and Frontier Research Project of Chongqing under Grant No. cstc2018jcyjAX0144, the Science and Technology Research Program of Chongqing Municipal Education Commission under Grant No. KJZD-M202000502, the Program of Chongqing Graduate Research and Innovation Project under Grant CYS19290, the Open Project of Key Laboratory No.CSSXKFKTM202007, Mathematical College, Chongqing Normal University, the Doctor Start-up Funding of Chongqing University of Posts and Telecommunications under Grant A2016-80. The authors thank all anonymous reviewers for their helpful comments and valuable suggestions.

Conflict of interest

The authors declare there is no conflict of interest.