Optimal harvesting strategy of a stochastic $n$-species marine food chain model driven by Lévy noises

1.   Introduction

Overfishing occurs when fish are harvested at a rate faster than the resource can replenish itself. This practice has consequences that extend beyond the depletion of single-species fish populations, potentially leading to global food insecurity, ecological extinction, and other devastating effects^[1,2,3,4]. Marine ecosystems and their associated livelihoods worldwide are threatened by overfishing. Nicholas et al. ^[5] have shown that overfishing is pushing over one-third of all sharks and rays towards global extinction. They suggest that marine life could face a future similar to that of land biodiversity, where human pressures have resulted in the loss of numerous species, possibly triggering a sixth mass extinction. According to IUCN Red List criteria ^[6], they estimate that one-quarter of 1,041 chondrichthyan fishes are threatened due to overfishing (both targeted and incidental). Therefore, it is urgent to establish the optimal harvesting strategy (OHS) for overfishing to prevent extinctions and promote population recovery. Consequently, investigating OHS has significant implications for marine ecosystems.

Mathematical models are widely utilized to analyze the OHS of populations ^{[7,8,9,10,11,12,13]}. Clark ^[7] proposes several deterministic population systems with harvesting and establishes the OHS. Martin and Ruan ^[8] investigate the predator-prey models with delay and prey harvesting. However, nature is full of random phenomena, and various populations in the ecosystem experience different forms of random perturbations. May ^[9] indicates that environmental noise affects competition coefficient, environmental capacity, the growth rate and other parameters of population systems. Lande et al.^[10] demonstrate that the extinction of numerous species caused by over-harvesting may be largely due to the ignorance of randomness. Therefore, it is crucial to investigate the OHS of population system in random environment. Li and Wang ^[11] establish a general stochastic logistic population model with white noises and explore the optimal harvesting policy. Wang et al.^[12] propose a stochastic predator-prey model with time delays and shows how randomness influences the extinction and OHS. Therefore, exploring the OHS of stochastic marine population systems is meaningful.

Lévy noise, another type of environmental perturbations, is considered below. During the process of species evolution, population growth is often affected by random disturbances in the environment, such as toxic pollutants, epidemics, earthquakes, floods and other factors ^{[14,15,16,17]}. For example, the Sandoz chemical accident in 1986 resulted in the death of over 600,000 fish in the Rhine ^[14], while the 2010 Gulf of Mexico oil spill caused a devastating impact on marine life in the Gulf of Mexico ^[15]. Bao et al. ^[16] propose that non-Gaussian Lévy noise is suitable for concluding these phenomena and revealed that Lévy noise can inhibit the explosive increase of population. Zhang et al. ^[17] also proposed to use Lévy noise to describe the random perturbations.

A lot of work has been done to study the stochastic food chain population dynamics model so far. Our previous work has investigated the dynamics of a stochastic three species food chain model in ^[18], the dynamics of a four-species stochastic food chain model with delays in ^[19] and the dynamic behaviors of a stochastic $n$-species food chain model with Lévy noises and white noises in ^[20] as follows:

By introducing a new research method, that is algebraic approach, the sufficient and necessary criteria for the stochastic extinction, persistence in mean and global stability are established in ^[20]. These conclusions directly generalize some research results on three and four species stochastic food chain dynamics models in ^[18,19]. However, due to the high dimensionality of model (1.1) and technical difficulties, the stability in distribution and optimal harvesting problems have not been discussed until now. In this paper, we consider the effects of harvesting and Lévy noises, then address these issues by proposing a stochastic $n$-species food chain model as follows:

On the basis of the research results in ^[20], our primary objective is to establish criteria for the stability in distribution of model (1.2). Additionally, we aim to establish criteria for the optimal harvesting strategy in the sense of sustainable yield. Lastly, we will verify the theoretical results through numerical examples and propose some open questions and propositions.

This paper is organized as follows. In Section 2, we provide the basic assumption, and further introduce several important lemmas that will be used to prove the major results. In Section 3, we state and prove the stability in distribution of model (1.2). In Section 4, we outline the non-existence and existence criteria for OHS and MESY. Finally, in Section 5, we illustrate our major results through numerical examples.

2.   Preliminaries

Denote by $(\Omega, \{\mathcal{F}_t\}_{t\geq0}, P)$ the complete probability space, and its filtration $\{\mathcal{F}_{t}\}_{t\geq0}$ satisfies the usual conditions. In this paper, we assume that model (1.2) is defined on the probability space $(\Omega, \{\mathcal{F}_t\}_{t\geq0}, P)$. In addition, for an integrable function $f(t)$ defined for $t\in [0, +\infty)$, we define the mean value of $f(t)$ on interval $[0, t]$ for any $t > 0$ by $\langle f(t)\rangle = \frac{1}{t}\int_0^tf(s)ds$.

In models (1.1) and (1.2), parameter $r_1 > 0$ is the intrinsic growth rate of species $x_1$, $r_i\leq 0$ represents the death rate of species $x_i$ for $i = 2, 3, \cdots, n$, $h_i\geq 0$ represents the harvesting rate of human for the species $x_i$ for $i = 1, 2, \cdots, n$, $x_i(t^-)$ is the left limit of $x_i(t)$ at time $t$, $a_{ii} > 0\; (i = 1, 2, \cdots, n)$ is the density dependent coefficient of species $x_i$, $a_{ii+1}\geq0\; (i = 2, 3, \cdots, n-1)$ is the capture rate of which species $x_{i+1}$ preys on species $x_i$, $a_{ii-1}\geq 0\; (i = 2, 3, \cdots, n)$ stands for efficiency of food conversion from species $x_{i-1}$ to species $x_i$. $B_i(t)\; (i = 1, 2, \cdots, n)$ are independent standard Brownian motion, and $\sigma_i^2 > 0$ is the intensity of $B_i(t)$. Lévy noises $\widetilde{N}(dt, du) = N(dt, du)-\lambda(du)dt$, where $N(dt, du)$ is a poisson counting measure with characteristic measure $\lambda$ on a Borel-measurable subset $Z$ of $(0, \infty)$ satisfying $\lambda(Z) < +\infty$. Functions $\gamma_i(u): Z\to R\; (i = 1, 2, \cdots, n)$ represent the effects of random jumps for species $x_i$ that are assumed to be continuous on $Z$. In this paper we further assume that Lévy noises $\widetilde{N}(dt, du)$ and Brownian motion $B_i(t)\; (i = 1, 2, \cdots, n)$ are also independent.

For any solution $(x_1(t), \cdots, x_n(t))$ of model (1.2), the initial value is given as follows.

where $x_{i0}\; (i = 1, 2, \cdots, n)$ are positive constants.

Denote $b_i = r_i+\int_{z}[\ln(1+\gamma_i(u))-\gamma_i(u)]\lambda(du)-\frac{\sigma_i^2}{2}$ and $c_i = b_i-h_i$ for $i = 1, 2, \cdots, n$, and $\alpha_k = (c_1, c_2, \cdots, c_k)^T$ for $k = 1, 2, \cdots, n$. Furthermore, we denote

Particularly, we have

Denote by $H_k$ the determinant of $A_k$. Let $A_k^{(i)}$ be the matrix gained from matrix $A_k$ by changing its $i$-th column to $\alpha_k$, and $D_k^{(i)}$ is the determinant of $A_k^{(i)}$. Furthermore, for any square matrix $P$ we denote by $P(i, j)$ the algebraic complement corresponding to the element of $i$-th row and $j$-th column of $P$.

For the Lévy noises in models (1.1) and (1.2), we always assume that the following condition is satisfied.

$\bf(H)$ For $i = 1, 2, \cdots, n$ and $u\in Z$, $1+\gamma_i(u) > 0$ holds, and a constant $z > 0$ exists such that

Using the same method as in ^[20], the following conclusions can be established for model (1.2).

Lemma 1. $H_k > 0$ for every $k = 1, 2, \cdots, n$.

Lemma 2. For any initial values $(x_{10}, x_{20}, \cdots, x_{n0})\in R_+^n$, model (1.2) possesses a globally unique positive solution $x(t) = (x_1(t), \cdots, x_{n}(t))\in R^n_+\; a.s.$ for all $t\geq 0$ satisfying the initial condition (2.1). Furthermore, for any constant $q > 0$, there are constants $C_i(q) > 0\; (i = 1, 2, \cdots, n)$ such that

Similar to ^[20], we define

Lemma 3. For any positive solution $(x_1(t), \cdots, x_n(t))$ of model (1.2), assume that $D_n^{(n)} > 0$ and $(1+\frac{a_{nn-1}}{a_{nn}}N_{n-1}) > 0$ if $n\geq3$, then $\lim_{t\to\infty}\langle x_i(t)\rangle = \frac{D_{n}^{(i)}}{H_n}$ a.s., $i = 1, 2, \cdots, n$.

A simple calculation we can obtain that for any integers $k, i, j\in\{1, 2, \cdots, n\}$ with $i, j\leq k$

Therefore, we further have

3.   Stability in distribution

Theorem 1. Let $(x_1(t; \psi), x_2(t; \psi), \cdots, x_n(t; \psi))$ and $(z_1(t; \psi^*), z_2(t; \psi^*), \cdots, z_n(t; \psi^*))$ be two solutions of model (1.2) with initial values $\psi, \psi^*\in R_+^n.$ Assume that there are positive constants $w_i\; (i = 1, 2, \cdots, n)$ such that

Then we have

Proof. We only need to show

Define the functions as follows

Applying Itô formula, then

where $a_{10} = a_{nn+1} = 0$ are stupulated. Define the following function

From (3.1) we obtain

where $a_{01} = a_{n+1n} = 0$ and $\omega_0 = \omega_{n+1} = 0$ are stipulated. Therefore,

which implies

Define the functions

Then, for any $t_1, t_2\in[0, +\infty)$ and each $i = 1, 2, \cdots, n$,

Applying Itô formula to model (1.2), then

For any $t_2 > t_1$ and $q > 1$, from the first equation of (3.6) and Hölder's inequality, then

Similarly,

and

In view of Theorem 7.1 in ^[21] on the Burkholder-Davis-Gundy inequality, for $1 < q\leq2$ and any $t_2 > t_1$, then

Using the Kunita first inequality in ^[21], there exists a constant $K_q > 0$ and from assumption $(H)$ we can obtain

From Lemma 2, there exists $C_i^{**}(q) > 0$ for $i = 1, 2, \cdots, n$ such that $\sup_{t\geq{-\gamma}}E[x_i^q(t)]\leq C_i^{**}(q)$. Therefore, from (3.7)–(3.10) there is a $\delta > 0$ satisfies that for $1 < q\leq2$ and any $t_1\geq0, t_2\geq0$ with $|t_2-t_1|\leq\delta$, such that

where

Similarly, we also obtain

for any $t_1\geq0, t_2\geq0$ with $|t_2-t_1|\leq\delta$ and $1 < q\leq2$. Thus, from (3.5), we obtain

Using similar argument, for $F_i(t)\; (i = 2, 3, \cdots, n)$, we also can obtain that there is a $\delta > 0$ for any $t_1\geq0, t_2\geq0$ with $| t_2-t_1|\leq\delta$ and $1 < q\leq2$

where

From (3.12) and (3.13), we obtain that $F_i(t), \; i = 1, 2, \cdots, n$ for $t\in(0, \infty)$ are uniformly continuous. Therefore, from (3.4) and Barbalat lemma in ^[22] we finally obtain that (3.1) holds. This completes the proof. □

Let $\mathcal{P}(R_+^n)$ be the space of all probability measures on $C(R_+^n)$. For $Q_1, Q_2\in\mathcal{P}(R_+^n)$, we define

where

Denote by $p(t, \phi, dx)$ the transition probability of process $x(t) = (x_1(t), x_2(t), \cdots, x_n(t))$. We have the following results.

Theorem 2. Assume that there exists constants $w_i > 0\; (i = 1, 2, \cdots, n)$ such that

Then model (1.2) is asymptotically stable in distribution, that is, there is a probability measure $v(\cdot)$ satisfies that for any initial function $\psi\in R_+^n$, the transition probability $p(t, \psi, \cdot)$ of solution $(x_1(t, \psi), x_2(t, \psi), \cdots, x_n(t, \psi))$ such that

Proof. For any initial data $\psi\in R_+^n$, the family probability $p(t, \psi, \cdot)$ is tight which obtained from Theorem 1 and Chebyshev's inequality. For any $f\in LB$ and $t, s > 0$, we compute

It is clear that there exists a constant $T > 0$ from (3.1), such that for all $t\geq T$,

Consequently

Since $f$ is arbitrary, we gain

That is,

Therefore $\{p(t, \psi, \cdot):t\geq0\}$ is Cauchy in the space $\mathcal{P}(R_+^n)$. Hence, there is a unique $v(\cdot)\in\mathcal{P}(R_+^n)$ such that

where $\psi_0 = (c_{10}, c_{20}, ..., c_{n0})^T\in R_+^n$, and $c_{10}, \cdots, c_{n0}$ are positive constants. It then follows from (3.1) that

Hence,

This completes the proof. □

4.   Optimal harvesting strategy

Since $H_n > 0$ in Lemma 1, the inverse matrix $A_n^{-1}$ exists. Let $B = (b_1, b_2, \cdots, b_n)^T$ and $\mathcal{H} = (h_1, h_2, \cdots, h_n)^T$. Moreover, let $E$ be the unit matrix, and $\mathcal{H}^* = (A_n(A_n^{-1})^T+E)^{-1}B\triangleq(h_1^*, h_2^*, \cdots, h_n^*)^T$.

Theorem 3. Assume that the following conditions hold

(1) $(1+\frac{a_{nn-1}}{a_{nn}}N_{n-1}) > 0$ if $n\geq 3$;

(2) there exists constants $w_i > 0\; (i = 1, 2, \cdots, n)$ such that

Then the following conclusions hold.

($\mathcal{A}_1$) If $D_n^{(n)}|_{h_i = h_i^*} > 0$, $A_n^{-1}+(A_n^{-1})^T$ is positive semi-definite and $h_i^*\geq0$ $(i = 1, 2, \cdots, n)$, then model (1.2) has the OHS $H = H^*$, and

($\mathcal A_2$) If any of the following conditions stands, then the OHS for model (1.2) does not exist.

$(\mathcal B_1)$ $D_n^{(n)}|_{h_i = h_i^*, i = 1, 2, \cdots, n}\leq0$;

$(\mathcal B_2)$ There is a $h_i^*$ with $h_i^* < 0$;

$(\mathcal B_3)$ $A_n^{-1}+(A_n^{-1})^T$ is not positive semi-definite.

Proof. Let $\mathcal{Y} = \{\mathcal{H} = (h_1, h_2, \cdots, h_n)^T\in R^n: D_n^{(n)} > 0, h_i\geq0, i = 1, 2, \cdots, n\}$, which means that for any $\mathcal{H}\in\mathcal{Y}$, the conclusion of Lemma 3 stands.

First, we prove the conclusion $(\mathcal{A}_1)$. Since $\mathcal{H}^*\in\mathcal{Y}$, the set $\mathcal{Y}$ is not empty. A unique probability measure $v(\cdot)$ is gained from Theorem 2. Then, $v(\cdot)$ is also strong mixing by Corollary 3.4.3 ^[23]. From Theorem 3.2.6^[23], we gain that the measure $v(\cdot)$ is ergodic. For any initial value $\psi\in R_+^n$, model (1.2) possess a positive gloabl solution $x(t) = (x_1(t), x_2(t), \cdots, x_n(t))$. It yields by Theorem 3.3.1 ^[23] that

where $\mathcal{H} = (h_1, h_2, \cdots, h_n)^T\in\mathcal{Y}$. Let $\varrho(z)$ be the stationary probability density of model (1.2), then we have

In view of the invariant measure is sole, there also exists a one-to-one correspondence among $\varrho(z)$ and its corresponding invariant measure. Hence, we have

Combining (4.2)–(4.4) with the conclusion of Lemma 3, then

Moreover, $S(\mathcal{H}) = \mathcal{H}^T(A_n)^{-1}(B-\mathcal{H})$ is derived by careful calculation. Computing the gradient of $S(\mathcal{H})$, then

Since $\frac{\partial \mathcal{H}^T}{\partial\mathcal{H}} = E$, then

By solving $\frac{\partial S(\mathcal{H})}{\partial \mathcal{H}} = 0$, the following critical value is gained

That is,

Calculating the Hessian matrix of $S(\mathcal{H})$, we obtain

Since $(A_n)^{-1}+{((A_n)^{-1})}^T$ is positive semi-definit, from the existence principle of extremum value of multi-variable function we obtain that $S(\mathcal{H})$ has a global maximum value at $\mathcal{H} = \mathcal{H}^*$. Obviously, $\mathcal{H}^*$ is also sole. Therefore, if $\mathcal{H}^*\in\mathcal{Y}$, that is $h^*_i\geq0\; (i = 1, 2, \cdots, n)$ and $D_n^{(n)}|_{h_i = h^*_i, i = 1, 2, \cdots, n} > 0$, then we finally obtain that $\mathcal{H}^*$ is an OHS, and MESY shown in (20).

Next, we prove the conclusion $(\mathcal{A}_2)$. We first assume that $(\mathcal{B}_1)$ or $(\mathcal{B}_2)$ stands. Suppose that $\widetilde\Gamma = (\gamma_1, \gamma_2, \cdots, \gamma_n)$ is the OHS for model (1.2), thus $\Gamma\in\mathcal{Y}$, that is

On the other hand, if $\Gamma\in\mathcal{Y}$ is the OHS, then $\Gamma$ is a unique solution of equation $\frac{\partial S(\mathcal{H})}{\partial\mathcal{H}} = 0$. Therefore, we have $(h_1^*, h_2^*, \cdots, h_n^*) = (\gamma_1, \gamma_2, \cdots, \gamma_n)$. Thus, condition (4.6) becomes into

However, this is impossible.

Finally, consider $(\mathcal B_3)$. It can be assumed that $(\mathcal{B}_1)$ and $(\mathcal{B}_2)$ fail to stand in this case. Therefore, $h_i^*\geq0\; (i = 1, 2, \cdots, n)$ and $D_n^{(n)}|_{h_i = h_i^*, i = 1, 2, \cdots, n} > 0$. Thus, $\mathcal{Y}$ is not empty, i.e., $(4.2)$–$(4.4)$ hold. Denote $(A_n)^{-1}+((A_n)^{-1})^T = (h_{ij})_{n\times n}$. Then, by calculating we have

Obviously, $h_{11} > 0$, which means that $(A_n)^{-1}+((A_n)^{-1})^T$ is not negative semi-definite. It can be further derived from the conditions in $(\mathcal{B}_3)$ that $(A_n)^{-1}+((A_n)^{-1})^T$ is indefinite. Hence, there is not the OHS if $\mathcal{B}_3$ holds. This completes the proof. □

Remark 1. In Theorem 3 we see that the conditions (1) and (2) are required. From the following argumentations we will find that when $n = 3$ and $n = 4$ condition (2) can bring out condition (1). Therefore, we propose the following proposition. However, when $n > 4$ whether condition (2) can bring out condition (1) still is an open problem. Thus, the following conjecture can be proposed.

Proposition 1. When $n = 3$ and $n = 4$, if condition (2) holds then condition (1) also holds.

Proof. Let $k = 3$, then condition (2) becomes to

From (2.2), by calculating we have

and

From (4.7) we derive

Therefore, we further have

This shows that $1+\frac{a_{32}}{a_{33}}N_2 > 0$.

Let $k = 4$, then condition (2) becomes to

From (2.2), by calculating, we have

and

From (4.8) we obtain

Therefore,

Since

we further obtain

This shows that $1+\frac{a_{43}}{a_{44}}N_3 > 0$. □

Conjecture 1. When $n > 4$, if condition (2) holds then condition (1) also holds.

5.   Numerical examples

To support our analysis results, we will provide the numerical examples for model (1.2). The numerical simulations can be performed qualitatively, and parameters are selected in an epidemiologically feasible way.

Example 1. In model (1.2) we take $n = 4$, $r_1 = 1.5$, $r_2 = -0.6$, $r_3 = -0.08$, $r_4 = -0.01$, $\gamma_1 = 0.1$, $\gamma_2 = 0.15$, $\gamma_3 = 0.15$, $\gamma_4 = 0.1$, $a_{11} = 1.6$, $a_{12} = 0.2$, $a_{22} = 2.2$, $a_{21} = 2.5$, $a_{23} = 0.1$, $a_{32} = 2.5$, $a_{33} = 2$, $a_{34} = 0.3$, $a_{43} = 1.3$, $a_{44} = 2$, $\sigma_1 = 0.1$, $\sigma_2 = 0.2$, $\sigma_3 = 0.1$ and $\sigma_4 = 0.1$

Choose $w_1 = 2.5$, $w_2 = 1.3$, $w_3 = 0.8$ and $w_4 = 0.5$, by calculating we have $w_1a_{11}-w_2a_{21} = 0.7500 > 0$, $w_2a_{22}-w_1a_{12}-w_3a_{32} = 0.3600 > 0$, $w_3a_{33}-w_2a_{23}-w_4a_{43} = 0.8200 > 0$ and $m_4a_{44}-m_3a_{34} = 0.7600 > 0$. Proposition 1 shows that $1+\frac{a_{43}}{a_{44}}N_3 > 0$. Thus, all four species in model (1.2) are persistent in mean which is shown in Figure 1 below. Furthermore, $h_1^* = 0.3250$, $h_2^* = 0.3914$, $h_3^* = 0.3264$, $h_4^* = 0.1281$ and $D_{4}^{(4)}|_{h_i = h_i^*, i = 1, 2, 3, 4} = 0.5201 > 0$ are obtained by calculating, all leading principal minors of $A_4+(A_4^{-1})^T$ are positive, which shows that the $A_4+(A_4^{-1})^T$ is positive semi-definite. Thus, the condition of conclusion ($\mathcal{A}_1$) is satisfied, then the OHS: $H^* = (0.3250\text{, }0.3914\text{, }0.3264\text{, }0.1281)^T$ and MESY$= 0.4321$ are gained, as Figure 1 shows below.

Example 2. In model (1.2) we take $n = 4$, $r_1 = 1.5$, $r_2 = -0.5$, $r_3 = -0.03$, $r_4 = -0.01$; $\gamma_1 = 0.15$, $\gamma_2 = 0.2$, $\gamma_3 = 0.1$ and $\gamma_4 = 0.05$. The rest parameters in model (1.2) are shown in the following three cases.

Case 2.1. $a_{11} = 1.6$, $a_{12} = 0.2$, $a_{22} = 2$, $a_{21} = 2.5$, $a_{23} = 1$, $a_{32} = 2.5$, $a_{33} = 2.1$, $a_{34} = 0.2$, $a_{43} = 0.1$, $a_{44} = 2$, $\sigma_1 = 0.1$, $\sigma_2 = 0.2$, $\sigma_3 = 0.2$ and $\sigma_4 = 0.1$.

Case 2.2. $a_{11} = 1.8$, $a_{12} = 0.2$, $a_{22} = 2$, $a_{21} = 1.5$, $a_{23} = 0.2$, $a_{32} = 1.5$, $a_{33} = 2.1$, $a_{34} = 0.8$, $a_{43} = 1.1$, $a_{44} = 2.5$, $\sigma_1 = 0.2$, $\sigma_2 = 0.2$, $\sigma_3 = 0.3$ and $\sigma_4 = 0.1$.

Case 2.3. $a_{11} = 0.8$, $a_{12} = 0.2$, $a_{22} = 2$, $a_{21} = 1$, $a_{23} = 0.15$, $a_{32} = 1.5$, $a_{33} = 2$, $a_{34} = 0.1$, $a_{43} = 1$, $a_{44} = 0.08$, $\sigma_1 = 0.2$, $\sigma_2 = 0.3$, $\sigma_3 = 0.3$ and $\sigma_4 = 0.1$.

For Case 2.1, choose $w_1 = 2.5$, $w_2 = 1.3$, $w_3 = 0.8$ and $w_4 = 1.4$, then we have $w_1a_{11}-w_2a_{21} = 0.75 > 0$, $w_2a_{22}-w_1a_{12}-w_3a_{32} = 0.10 > 0$, $w_3a_{33}-w_2a_{23}-w_4a_{43} = 0.24 > 0$ ºÍ $m_4a_{44}-m_3a_{34} = 2.64 > 0$. Furthermore, $h_1^* = 0.3621$, $h_2^* = 0.4773$, $h_3^* = 0.4665$, $h_4^* = 0.0051$ are gained by calculating. Meanwhile, $D_{4}^{(4)}|_{h_i = h_i^*, i = 1, 2, 3, 4} = -0.1475 < 0$. Thus the condition ($\mathcal{B}_1$) of ($\mathcal{A}_2$) is satisfied. Therefore, the OHS do not exist. The numerical simulations are given in Figure 2.

For Case 2.2, choose $w_1 = 2.2$, $w_2 = 1.1$, $w_3 = 1.1$ and $w_4 = 1.4$, then we have $w_1a_{11}-w_2a_{21} = 2.31 > 0$, $w_2a_{22}-w_1a_{12}-w_3a_{32} = 0.11 > 0$, $w_3a_{33}-w_2a_{23}-w_4a_{43} = 0.55 > 0$ ºÍ $m_4a_{44}-m_3a_{34} = 2.62 > 0$. Furthermore, $h_1^* = 0.9445 > 0$, $h_2^* = 0.2174 > 0$, $h_3^* = 0.0039 > 0$ and $h_4^* = -0.0124 < 0$ are obtained by calculating. Thus the condition ($\mathcal{B}_2$) of conclusion ($\mathcal{A}_2$) is satisfied. Therefore, the OHS do not exist. The numerical simulations are given in Figure 3.

For Case 2.3, choose $w_1 = 2.8$, $w_2 = 2$, $w_3 = 2$ and $w_4 = 2.8$, then we have $w_1a_{11}-w_2a_{21} = 0.24 > 0$, $w_2a_{22}-w_1a_{12}-w_3a_{32} = 0.44 > 0$, $w_3a_{33}-w_2a_{23}-w_4a_{43} = 0.90 > 0$ ºÍ $m_4a_{44}-m_3a_{34} = 0.024 > 0$. Furthermore, by calculating we obtain that the fourth leading principle minor of $A_4+(A_4^{-1})^T$ is equal to $-6.2127 < 0$, which means that $A_4+(A_4^{-1})^T$ is negative semi-definite. Thus the condition ($\mathcal{B}_3$) is satisfied. Therefore, the optimal harvesting strategy do not exist. The numerical simulations are given in Figure 4.

6.   Conclusions

In this paper, based on the effect of harvesting and Lévy noises, a stochastic marine $n$-species system is investigated. On the basis of the research results in ^[23], we first establish the criterion for the stability in distribution of model (1.2). Then, the criteria for existence and nonexistence of OHS and the corresponding MESY are derived. This shows that the environmental noises affect the OHS and MESY.

Compared with the existing results, our first contribution is to investigate the stochastic stability and optimal harvesting problems for a stochastic marine $n$-species food chain model with Lévy noises and white noises. Second, the criteria on the stability in distribution and optimal harvesting strategy for stochastic $n$-species food chain model are established for the first time. Third, the noises intensity can easily affect the dynamics and MESY of marine populations, and further lead to the imbalances of marine ecology. Finally, we find that for better manage and protect marine resources, the establishment of an optimal harvesting strategy should fully consider the impact of noise intensity.

There are still some remained problems that need to be followed up. For example, in Conjecture 1, we assume if $w_1a_{11}-w_2a_{21} > 0\text{, } \; w_ia_{ii}-w_{i-1}a_{i-1i}-w_{i+1}a_{i+1i} > 0\; (i = 2\text{, }3\text{, }\cdots\text{, } n-1)\text{, }\; w_na_{nn}-w_{n-1}a_{n-1n} > 0$, then $1+\frac{a_{nn-1}}{a_{nn}}N_{n-1} > 0$, how to prove is an open problem. Second, marine population are easily affected by climate change, consider marine population model with regime switching is necessary. Finally, more complex and realistic stochastic system with nonlinear functional response, delays and Allee effect need to be investigated. These issues may be investigated in the future.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (Grant No. 12061079), Foundation of Fujian Educational Committee (Grant No. JAT220189), Startup Foundation of Jimei University(Grant No.ZQ2022025).

Conflict of interest

The authors declare there is no conflicts of interest.