Dynamic analysis and validation of a prey-predator model based on fish harvesting and discontinuous prey refuge effect in uncertain environments

1.   Background and motivation

As an important branch of ecosystem research, the study on population dynamic model has always been one of the important topics in bioscience and mathematics. Especially in the fishery industry, understanding the relations between predator fish and prey fish is helpful for the utilization and development of fish resources in a sustainable and reasonable way. On the theoretical side, mathematical models play a key role in understanding changes in biological systems and the effects of related control measures ^[1,2]. The classical model in literature describing the relationship between predator species and prey species is the well-known Lotka-Volterra model ^[3,4]. Subsequently, for different application scenarios, scholars considered different factors in the modeling process ^[5,6]. In natural ecosystems, most prey species perceive the danger of predators well and hide in cover to avoid being eaten when predators are present. Thus, the predator does not use the entire prey population as its food resource, and the prey refuge concept was introduced into the predator-prey system ^[7]. In the literature, various types of prey refuge effect were studied, which can be roughly divided into three types: constant quantity ^[7], constant proportion ^[8], and variable proportion ^[9]. In this paper, inspired by the piecewise form of prey refuge, a Filippov type prey-predator model with piecewise form of prey refuge is put forward, where prey refuge makes sense once the number of predator exceeds a certain threshold related to the number of prey. The Filippov theory, a primary tool for analyzing Filippov models, is widely applied and studied. For example, Tang and Liang ^[10] investigated a Filippov predator-prey system with non-smooth refuge. Chen and Huang ^[11] analyzed a Filippov ratio-dependent prey-predator model with behavioral refuges caused by prey instinct anti-predator behavior. Li et al. ^[12] analyzed a Filippov-type plant disease models with an interaction ratio threshold. Subsequently, Li et al. ^[13] studied a Filippov predator-prey model with two thresholds for integrated pest management. In this study, considering the piecewise form of prey refuge, a fishery model involving discontinuous prey refuge effect is presented and investigated.

In the actual world, the natural environment is always in change and fluctuation, which will naturally influence the survival and reproduction of the species living there. To describe the influence caused by environment fluctuation, scholars introduced different models with interval-valued imprecise parameters and fuzzy parameters, for example, Pal et al. ^[14] replaced the deterministic parameters in the predator-prey models with interval-valued imprecise parameters and fuzzy parameters. Zhang and Zhao ^[15] discussed the bifurcation and optimal harvesting of a diffusive predator-prey system with delays and interval biological parameters. Pal et al. ^[16,17] discussed the stability and bionomic analysis of a prey-predator with fuzzy parameters. Xiao et al. ^[18] analyzed a competition fishery model with interval-valued parameters and discussed the extinction, coexistence, bionomic equilibria, and optimal harvesting policy. Wang et al. ^[19] incorporated prey refuge into a predator-prey system with imprecise parameter estimates. Yu et al. ^[20] analyzed a predator-prey fishery model with interval imprecise parameters, taking into account the interaction between predators in the system and the refuge effect of prey. Meng and Wu ^[21] analyzed the dynamics of a fuzzy phytoplankton-zooplankton model with refuge, fishery protection, and harvesting. Wang et al. ^[22] discussed the stability and optimal harvesting of a predator-prey system combining prey refuge with fuzzy biological parameters. Chen and Zheng ^[23] discussed the diffusion-driven instability of a predator-prey model with interval biological coefficients. Xu et al. ^[24] discussed the optimal harvesting of a fuzzy water hyacinth-fish model with Kuznets curve effect. Guo et al. ^[25] analyzed the dynamics of two fishery capture models with a variable search rate and fuzzy biological parameters. Cao et al. ^[26] discussed the Hopf bifurcation in a predator-prey model under fuzzy parameters involving prey refuge and fear effects. Studies on freshwater fish have shown that fish species are more vulnerable to environmental changes. In addition, environmental changes affect different fish species in different ways ^[27]. In view of this feature, species related imprecise parameters were introduced into the prey-predator model ^[28] and the impact of imprecise parameters on the dynamics of the systems were discussed. In the current work, a Filippov type prey-predator model with triangle fuzzy parameters is studied and the impact of fuzzy parameters on the dynamic behaviour is analyzed.

Fishing is the main way for human to obtain fish resources, and its pattern presents two types: continuous and discontinuous. For continuous fishing activities, the modeling is relatively simple, which can be directly described by adding fishing items in the model. For intermittent fishing patterns, fishing activities are usually taken at discrete moments. Among many intermittent fishing activities, state-dependent feedback fishing is a typical one, which takes into account the current state of prey or predator and is able to avoid destroying the sustainability of fish resources. There are many studies on state-dependent feedback control in the literature, which can be roughly divided into several types, such as prey-dependent ^[29], predator-dependent ^[30], nonlinear-dependent ^[31], weighted-dependent ^[32], and ratio-dependent ^[33]. Among these kinds of state-dependent feedback control, ratio-dependent feedback control considers the relation between prey population and predator population and implements control when the ratio of predator population to prey population reaches a certain threshold. In this study, a state's linearly dependent feedback fishing strategy is considered and the fishery model with such control is analyzed.

Motivated by above discussions, a Filippov type prey-predator model with discontinuous refuge effect, triangle fuzzy biological parameters and continuous harvesting is put forward and analyzed. The paper is organized as follows: In Section 2, two types of prey-predator models with different fishing model and fuzzy parameters are presented and followed by a presentation of some preliminaries. In Section 3, the dynamics of Filippov type prey-predator with a continuous harvesting strategy is analyzed. Then, the complex dynamics of the fishery model with state linearly dependent feedback harvesting is investigated. In Section 4, numerical simulations are presented to illustrate the main results. In the last section, a conclusion is summarized with a presentation of the future work.

2.   Model formulation and preliminaries

Prey refuge is a common phenomenon among fish species. Let $x_r$ denote the volume of prey in refuge. In this work, a discontinuous refuge effect is considered, that is,

where $m\in]0, 1[$ is the proportion of refuge prey; $y_T > 0$ is the minimum level of predators, when the density of predators is lower than $y_T$, prey will come out of the refuge; $n > 0$ is the threshold for the ratio of predator to prey in the system, and if the ratio exceeds the threshold of $n$ in the case of sufficient prey density, the prey will choose to hide in the shelter.

2.1.   Fishery model with continuous harvesting

Then, the Filippov type prey-predator model with continuous harvesting is expressed as

where $t\in [0, +\infty[$, $x(t)$ and $y(t)$ represent the biomass of prey and predator species at time $t$, respectively. $r$ characterizes the prey's growth rate; $K$ characterizes the prey's environmental capacity; $a$ characterizes the predator's predation rate, $x(t)-x_r$ is the density of prey outside the shelter; $d$ characterizes the predator's death rate; $c$ characterizes the conversion efficiency from prey biomass into predator biomass; $E_1$ and $E_2$ represent fishing effort for prey and predator, respectively; $q_1$ and $q_2$ represent the capture rate for prey and predator species.

To consider the impact of environmental changes and fluctuations, triangular fuzzy numbers (TFNs) ^[20] are adopted to describe the uncertainty of parameters. For a TFN $\tilde{U}\equiv(u_1, u_2, u_3)$ and $\alpha\in (0, 1]$, the $\alpha$-cut set is denoted by $[U_{l(\alpha)}, U_{r(\alpha)}]$, where $U_{l(\alpha)} = \{x:\mu_{\tilde{U}}(x)\geq\alpha\} = u_1+\alpha(u_2-u_1)$, $U_{r(\alpha)} = \{x:\mu_{\tilde{U}}(x)\geq\alpha\} = u_{3}+\alpha(u_{3}-u_{2})$. Given that the birth rate of prey, the death rate and conversion rate of predator most susceptible to environmental changes, these three parameters are assumed to present some imprecision, represented by TFNs, that is, $\tilde{r} = (r_L, r_M, r_R)$, $\tilde{d} = (d_L, d_M, d_R)$ and $\tilde{c} = (c_L, c_M, c_R)$. Using theory of $\alpha$-cut fuzzy number, we introduce fuzzy triangle parameters into above model:

Using the utility function method ^[16], there is

where $0\leq w_1, w_2\leq1$. For convenience, define

Combining Eqs (2.4) and (2.1) gives

For $y(t) < nx(t)+y_{T}$, system complies with the model

and for $y > nx+y_T$, system complies with the model

2.2.   Fishery model with threshold harvesting

In this subsection, we consider a scenario where fishing is allowed only when prey and predator populations exceed a certain limit. Let $y_H$ be the minimum fishing level of the predator, below which fishing activity may cause the extinction of the predator fish. Moreover, predator fish feed on prey populations, and when the ratio of predator to prey exceeds a certain threshold, denoted by $l$, fishing activity is conducive to the sustainable development of fish resources. Based on the above consideration, we establish the following harvesting model with uncertain parameters:

where $l > 0$, $y_H > 0$ are predetermined constants.

2.3.   Preliminaries

2.3.1.   TFN

Let $\tilde{U}$ be a fuzzy set on the real set $\mathbb{R}$, i.e., $\tilde{U}\in \mathcal{F}(\mathbb{R})$, $\mu_{\tilde{U}}(\cdot)$ be the membership function of $\tilde{U}$.

Definition 1 ($\alpha$-cut set ^[16,25]). For $\alpha \in ]0, 1]$, the $\alpha$-cut set for $\tilde{U}$ is defined as $\tilde{U}_{\alpha} = \left\{x: \mu_{\widetilde{U}}(x) \geq \alpha\right\}$.

Definition 2 (TFN^[16,25]). If $\tilde{U}$ is normal (i.e., there is $x\in \mathbb{R}$ and $\mu_{\tilde{U}}(x) = 1$), and for any $\alpha\in ]0, 1[$, $\tilde{U}_{\alpha}$ is a closed interval, then $\tilde{U}$ is said to be a fuzzy number (FN). A TFN $\tilde{U}\equiv (u_L, u_M, u_R)$ is a FN with membership defined by

Clearly, $\alpha$-cut set of TFN $\tilde{U}\equiv (a_L, a_M, a_R)$ is $[U_{l(\alpha)}, U_{r(\alpha)}]$, where $U_{l(\alpha)} = \inf\left\{x: \mu_{\tilde{U}}(x)\geq \alpha \right\} = u_L +\alpha (u_M - u_L)$ and $U_{r(\alpha)} = \sup\left\{x: \mu_{\tilde{U}}(u)\geq \alpha \right\} = u_R +\alpha (u_R - u_M)$.

Definition 3 (Utility function ^[16,25]). Given $U_i$, $i = 1, 2, \cdots, N$, and denoted $w_i$ as the weight of items $U_i$, $\Sigma^N_{i = 1} w_i = 1$. Then, a utility function $U$ is defined by $U = \Sigma^{N}_{i = 1} w_i U_i$.

2.3.2.   Filippov system

Consider a piecewise-continuous system

where $H: \mathbb{R}^2_{+}\rightarrow \mathbb{R}$, and the discontinuous demarcation is $\Sigma = \left\{(u, v) \in \mathbb{R}^2_{+}: H(u, v) = 0\right\}$.

Let

Then, $\mathbf{F}_{i}^{m} H = \langle\nabla\left(\mathbf{F}_{i}^{m-1} H\right), \mathbf{F}_{i}\rangle$ for $i = 1, 2$, $m \in \mathbb{N}$ with $m\geq 2$. The discontinuous demarcation $\Sigma$ can be distinguished into three regions:

1) sliding region: $\Sigma_{s} = \left\{(u, v) \in \Sigma: \mathbf{F}_{1} H < 0\right. \mbox{and} \left.\mathbf{F}_{2} H > 0\right\}$;

2) crossing region: $\Sigma_{c} = \left\{(u, v) \in \Sigma: \mathbf{F}_{1} H \cdot \mathbf{F}_{2} H > 0\right\}$;

3) escaping region: $\Sigma_{e} = \left\{(u, v) \in \Sigma: \mathbf{F}_{1} H > 0\right. \mbox{and} \left.\mathbf{F}_{2} H < 0\right\}$.

The dynamics of system (2.9) along $\Sigma_{s}$ is determined by

where $\mathbf{F}_{s} = \lambda \mathbf{F}_{1}+(1-\lambda) \mathbf{F}_{2}$ with $\lambda = \frac{\mathbf{F}_{2} H}{\mathbf{F}_{2} H-\mathbf{F}_{1} H} \in]0, 1[$.

Definition 4 (Real, virtual and pseudo-equilibrium ^[11]). For system (2.9), $E^{*}$ is a real equilibrium if $\exists i \in\{1, 2\}$ so that $\mathbf{F}_{i}\left(E^{*}\right) = 0$, $E^{*} \in S_{i}$; $E^{*}$ is a virtual equilibrium if $\exists i, j \in\{1, 2\}, i \neq j$ so that $\mathbf{F}_{i}\left(E^{*}\right) = 0$, $E^{*} \in S_{j}$; $E^{*}$ is a pseudo-equilibrium if $\mathbf{F}_{s}\left(E^{*}\right) = \lambda \mathbf{F}_{1}\left(E^{*}\right)+(1-\lambda) \mathbf{F}_{2}\left(E^{*}\right) = 0, H\left(E^{*}\right) = 0$, $\lambda = \frac{\mathbf{F}_{2} H}{\mathbf{F}_{2} H-\mathbf{F}_{1} H} \in]0, 1[$.

2.4.   Impulsive semi-continuous system

For a given planar model

Definition 5 (Order-$k$ periodic solution ^[32]). The solution $\tilde{\mathbf{z}}(t) = (\tilde{u}(t), \tilde{v}(t))$ of system (2.10) is called periodic if there exists $m (\geqslant 1)$ satisfying $\tilde{\mathbf{z}}_{m} = \tilde{\mathbf{z}}_{0}$. Furthermore, $\tilde{\mathbf{z}}$ is an order-$k$ T-periodic solution with $k\triangleq\min\{l|1\leq l\leq m, \tilde{\mathbf{z}}_{l} = \tilde{\mathbf{z}}_{0}\}$.

Lemma 1 (Analogue of Poincaré criterion ^[32]). The order-$k$ $T$-periodic solution $\mathbf{z}(t) = (\xi(t), \eta(t))^T$ of system (2.10) is orbitally asymptotically stable if $\vert \mu _k \vert < 1$, where

with

$f_1^+ = f_1(\xi(\tau_j^ +), \eta(\tau_j^ +))$, $f_2^+ = f_2(\xi(\tau_j^+), \eta(\tau_j^+))$ and $f_1$, $f_2$, $\frac{\partial {I_1} }{\partial u}$, $\frac{\partial {I_1} }{\partial v}$, $\frac{\partial {I_2} }{\partial u}$, $\frac{\partial {I_2} }{\partial v}$, $\frac{\partial \phi }{\partial u}$, $\frac{\partial \phi }{\partial v}$ are calculated at $(\xi(\tau_j), \eta(\tau_j))$, where $0 < \tau_1 < \cdots < \tau_{j-1} < \tau_j < \tau_{j+1} < \cdots < \tau_k = T$ and $\phi(\xi(\tau_j), \eta(\tau_j)) = 0$.

3.   Dynamic properties of systems (2.5)

Define

Obviously, when $E_1\geq \overline{E}_1$, there is $\frac{\text{d}x}{\text{d}t} < 0$, then $x(t)\rightarrow 0$ when $t\rightarrow \infty$.

To avoid the extinction of prey and predator fish populations by fishing activity}}, it is required that $E_1 < \overline{E}_1$ and $E_2 < \overline{E}_2(E_1)$.

3.1.   Positivity and bounded-ness of the solutions

Theorem 1. For given $w_1$, $w_2$, and $\alpha$, the solution of the Filippov system (2.5) with $x(0) = x_0 > 0$ and $y(0) = y_0 > 0$ always keeps positive, that is, $x(t) > 0$, $y(t) > 0$ for $t\in[0, +\infty)$.

Proof. Suppose that $(x(0), y(0))\in S_1 = \{(x, y)\in \mathbb{R}^2_{+}|y-nx-y_T < 0\}$. If $\mathbf{z}(t) = (x(t), y(t))$ intersects with region $\Sigma$ and then stay in region $\Sigma$ all the time, it can be easily obtained that $x(t) > 0$ and $y(t) > 0$ for $t\in[0, +\infty)$. Otherwise, define $t_i$, $i = 1, 2, ...$ as the time when trajectory intersects with $\Sigma$ and subsequently enters into another region. Note that due to the definition of function $x_r$, we can write

Since in $S_i$ there is

then we have

thus it has $x(t) > 0$ for $t\in[0, +\infty[$. Similarly, there is $y(t) > 0$ for $t\in[0, +\infty[$. □

Theorem 2. For given $w_1$, $w_2$, and $\alpha$, the solution $(x(t), y(t))$ of the Filippov system (2.5) with $x(0) = x_0 > 0$ and $y(0) = y_0 > 0$ is uniformly bounded.

Proof. Define $z(t) = \hat{c}x(t)+y(t)$. Then,

Choosing $0 < s\le \hat{d}+q_2E_2$, we have

which implies that

For $t\to +\infty$, there is $0\le z(t)\le \frac{\Theta}{s}$. Moreover, if $z_0 = z(x(0), y(0)) < \frac{\Theta}{s}$, then $0\le z(t)\le \frac{\Theta}{s}$ for all $t\geq 0$. □

3.2.   Dynamics of the subsystems (2.6) and (2.7)

Denote

Note that $P_0(0, 0)$ is a real equilibrium for subsystem (2.6) and a virtual equilibrium for subsystem (2.6). If $E_1 < \overline{E}_1$, then $P_{B_1}\left(\mathrm{K}^{1}_{\mathbf{E}}, 0\right)$ exists, which is a real boundary equilibrium for subsystem (2.6); $P_{B_2}\left(\mathrm{K}^{2}_{\mathbf{E}}, 0\right)$ exists, which is a virtual equilibrium for subsystem (2.7).

Theorem 3. For subsystems (2.6) and (2.7) with given $w_1$, $w_2$, and $\alpha$, $P_0(0, 0)$, $P_{B_1}\left(\mathrm{K}^{1}_{\mathbf{E}}, 0\right)$, and $P_{B_2}\left(\mathrm{K}^{2}_{\mathbf{E}}, 0\right)$ are unstable if $E_1 < \overline{E}_1$, $E_2 < \overline{E}_2(E_1)$. Moreover, subsystem (2.6) has a unique globally asymptotically stable interior equilibrium $P_1(x^*_1, y^*_1)$, subsystem (2.7) has a unique globally asymptotically stable interior equilibrium $P_2(x^*_2, y^*_2)$, where

Proof. In case of $E_1 < \overline{E}_1$, $E_2 < \overline{E}_2(E_1)$, the equation set $\left\{\begin{array}{c} f_{11}(x, y) = 0 \\ f_{12}(x, y) = 0 \end{array}\right.$ has a positive solution $(x, y) = (x^*_1, y^*_1)$, that is, subsystems (2.6) has a unique interior equilibrium $P_1(x^*_1, y^*_1)$. Similarly, the subsystem (2.7) has a unique interior equilibrium $P_2(x^*_2, y^*_2)$.

The Jacobian matrix of the subsystems (2.6) at $\bar{P}(\bar{x}, \bar{y})$ is

At $P_{0}(0, 0)$, there is

then $\lambda_1 = \hat{r}-q_1E_{1} > 0$ and $\lambda_2 = -d-q_2E_{2} < 0$, which implies that $O(0, 0)$ is saddle and unstable.

At $P_{B_1}\left(\frac{K(\hat{r}-q_1E_1)}{\hat{r}}, 0\right)$, there is

then $\lambda_1 = -\hat{r} < 0$ and $\lambda_2 = a\hat{c}\frac{K(\hat{r}-q_1E_1)}{\hat{r}}-d-q_2E_{2} > 0$, which implies that $P_{B_1}\left(\frac{K(\hat{r}-q_1E_1)}{\hat{r}}, 0\right)$ is saddle and unstable.

At $P^*_1(x^*_1, y^*_1)$, there is

then $\lambda_1+\lambda_2 = -\frac{\hat{r}}{K}x^*_{1}$, $\lambda_1\lambda_2 = a^2\hat{c}x^*_{1} y^*_{1} > 0$, which implies that $P^*_{1}(x^*_{1}, y^*_{1})$ is locally asymptotically stable.

Define $D(x, y) = \frac{1}{xy}$. Then, we have

then by the Bendixson-Dulac theorem ^[34], there does not exist a closed orbit in $\mathbb{R}^2_{+}$, so $P^*_{1}(x^*_{1}, y^*_{1})$ is globally asymptotically stable.

The stability of the equilibria $P_{0}$, $P_{B_2}$, and $P^*_2$ for the subsystems (2.7) can be proved similarly. □

3.3.   Sliding mode dynamics

It is obvious that $x^*_1 < x^*_2$ and $y^*_1 < y^*_2$. Define

Then, the slope and the intercept are, respectively,

Theorem 4. For Filippov system (2.5) with given $w_1$, $w_2$ and $\alpha$, for case-I): $y_T > y_l$, $n > (y_2-y_T)x_2$, $P_1$ is real, $P_2$ is virtual; for case-II): $y_T > y_l$, $0 < n < (y_1-y_T)x_1$, $P^*_1$ is virtual, $P^*_2$ is real; for case-III): $y_T > y_l$, $(y_1-y_T)x_1 < n < (y_2-y_T)x_2$, $P^*_1$ is real, $P^*_2$ is real; for case-IV): $y_T < y_l$, $n < (y_2-y_T)x_2$, $P^*_1$ is virtual, $P^*_2$ is real; for case-V): $y_T < y_l$, $n > (y_1-y_T)x_1$, $P^*_1$ is real, $P^*_2$ is virtual; for case-VI): $y_T < y_l$, $(y_2-y_T)x_2 < n < (y_1-y_T)x_1$, $P^*_1$ is virtual, $P^*_2$ is virtual.

Proof. For case-Ⅰ, that is, $y_T > y_l$, $n > \frac{y_2-y_T}{x_2}$, both of $P^*_1$ and $P^*_2$ are under the line $l_H$, as presented in Figure 1. $\mathbf{F}_1 (P^*_1) = 0$ and $P^*_1\in S_1$, so $P^*_1$ is a real equilibrium; $\mathbf{F}_2 (P^*_2) = 0$ and $P^*_2\in S_1$, so $P^*_2$ is a virtual equilibrium. Similarly, the results for case Ⅱ-Ⅵ can be proved.

□

Next, it discusses the existence of pseudo- equilibrium. Since $\nabla H = (n, -1)$ and

Since

then there exists a unique $x^1_s\in]0, \mathrm{K}_{E}[$ such that $M(x^1_s) = 0$, where

with

For $x\in ]0, x^1_s[$, there is $\mathbf{F}_1 H > 0$}}, and for $x\in ]x^1_s, K_E[$, there is $\mathbf{F}_1 H < 0$.

Similarly, there is

Since

then there exists a unique $x^2_{s}\in]0, \mathrm{K}^{1}_E[$ such that $G(x^2_s) = 0$, where

with

For $x\in]0, x^2_{s}[$, there is $\mathbf{F}_{2}H > 0$, and for $x\in]x^2_{s}, \mathrm{K}^{2}_{E}[$, there is $\mathbf{F}_{2}H < 0$.

Since $M(0) = G(0) = \hat{d}y_{T}+q_{2}E_{2}y_{T}$ and $M(x) < G(x)$ for $x > 0$, then it has $x_{s1} < x_{s2}$. Thus, $\sum_s = \{(x, y)\in\sum: x_{s1} < x < x_{s2}\}$ and $\sum_c = \sum_{c_1}\cup\sum_{c_2}$, where $\sum_{c_{1}} = \{(x, y)\in\sum\colon0 < x < x^*_{1}\}$ and $\sum_{c_{2}} = \{(x, y)\in\sum\colon x > x^*_{2}\}$.

Utilizing the Filippov convex method, the dynamics on $\sum_s$ is determined by

where

Since on $\sum_s$, there is $H(x, y) = 0$, then we have

Define

Then, we have the following result:

Theorem 5. For given $w_1$, $w_2$, and $\alpha$, if $\Delta > 0$ and $Q(x^1_s)Q(x^2_s) < 0$, then Filippov system (2.5) has a unique pseudo-equilibrium $E_P(x_P, nx_P+y_T)$.

3.4.   Complex dynamics system (2.8)

For system (2.8), we have

To ensure that $\mathcal{M}$ and $\mathcal{N}$ do not intersect, it requires that $q_2E_2 > q_1E_1$, in such case $\mathcal{N}$ is below $\mathcal{M}$. Let $L_1, L_2, L_3, L_4$ be the isolines in system, that is,

According to the relative position between $\mathcal{M}$, $\mathcal{N}$, and $l_H$, two situations are discussed for the dynamic of system (2.8):

Case Ⅰ: $n > l$, $y_T > y_H$, that is, both $\mathcal{M}$ and $\mathcal{N}$ lie below $l_H$ in $S_1$.

Case Ⅱ: $n < l(1-q_2E_2)/(1-q_1E_1)$, $y_T < y_H(1-q_2E_2)$, that is, both $\mathcal{M}$ and $\mathcal{N}$ lies above $l_H$ in $S_2$.

Definition 6 (Successor function). For a point $S\in \mathcal{N}$, if the trajectory from $S$ directly intersects $\mathcal{M}$, denote the intersection point by $S^{-}\in \mathcal{M}$. And then $S^{-}$ is impulsed to $S^{+}\in \mathcal{N}$ due to impulse effect. In such case, we can define $f^{\text{I}}_{\text{sor}}$: $\mathcal{N}\rightarrow \mathbb{R}$, $S\rightarrow f^{\text{I}}_{\text{sor}}(S)\triangleq y_{S^{+}}-y_{S}$. If the trajectory from $S$ first passes through $\mathcal{N}$, and then intersects $\mathcal{M}$, then we can define $f^{\text{II}}_{\text{sor}}$: $\mathcal{N}\rightarrow \mathbb{R}$, $S\rightarrow f^{\text{I}}_{\text{sor}}(S)\triangleq y_{S^{+}}-y_{S}$.

Theorem 6. For given $w_1$, $w_2$ and $\alpha$ and Case I, an order-1 periodic solution exists in system (2.8) for any one of the following conditions: I-1) $y_T > y_l$, $n > (y_2-y_T)/x_2$, $0 < y_H < H_1\triangleq y_1-lx_1$; I-2) $y_T > y_l$, $0 < n < (y_1-y_T)/x_1$; I-3) $y_T > y_l$, $(y_1-y_T)/x_1 < n < (y_2-y_T)/x_2$; I-4) $y_T < y_l$, $n < (y_2-y_T)/x_2$; I-5) $y_T < y_l$, $n > (y_1-y_T)/x_1$, $0 < y_H < H_1$; I-6) $y_T < y_l$, $(y_2-y_T)/x_2 < n < (y_1-y_T)/x_1$.

Proof. For case Ⅰ-1) $y_T > y_l, n > (y_2-y_T)/x^*_2$, $P^*_1$ is real equilibrium, $P^*_2$ is virtual equilibrium. Moreover, $P^*_1$ is locally asymptotically stable.

Denote $B$ as the interaction point between $L_1$ and $\mathcal{N}$. Since $0 < y_H < H_1\triangleq y_1-lx_1$, then $f_{\text{sor}}^{\text{I}}(B) < 0$. Let $A\in\mathcal{N}$ such that

The coordinates of $A$ are obtained from the following equations:

Define $\hat{E_2}\triangleq (1-y_A/y_{A^{-}})/q_2$. Then,

1) for $E_2 = \hat{E_2}$, there is $f_{\text{sor}}^{\text{I}}(A) = 0$, that is, the orbit from $A$ forms an order-1 periodic orbit;

2) for $E_2 < \hat{E_2}$, there is $f_{\text{sor}}^{\text{I}}(A) > 0$. Then, $\exists S'\in\overline{AB}\subset\mathcal{N}$ such that $f_{\text{sor}}^{\text{I}}(S') = 0$;

3) for $E_2 > \hat{E_2}$, there is $f_{\text{sor}}^{\text{I}}(A) < 0$. Then, $f_{\text{sor}}^{\text{II}}(A^{+}) > 0$. Besides, for $\varepsilon =$d$(A, A^+)/4 > 0$, $\exists \delta < \varepsilon$ and $A_1\in U(A, \delta)\cap\mathcal{N}$ so that $d(A^+, A_{1}^{+}) < \epsilon$, then it has $f_{\text{sor}}^{\text{II}}(A_1) < 0$. Thus $\exists S'\in\overline{A_1B^{+}}\subset\mathcal{N}$ such that $f_{\text{sor}}^{\text{I}}(S') < 0$.

To sum up, $\exists S'\in\mathcal{N}$, the trajectory of system (2.8) starting from $S'$ forms an order-l periodic solution. Similarly, the results for case Ⅰ-2)–Ⅰ-6) can be proved. □

Theorem 7. For given $w_1$, $w_2$, and $\alpha$ and Case II, an order-1 periodic solution exists in system (2.8) for any one of following two conditions: II-1) $y_T > y_l$, $0 < n < (y_1-y_T)/x_1$, $0 < y_H < H_2\triangleq y_2-lx_2$; II-2) $y_T < y_l$, $n < (y_2-y_T)/x_2$, $0 < y_H < H_2$.

Proof. The proof is similar to that of Theorem 5 and is omitted here. □

Let $\mathbf{z}(\mathrm{t}) = \left(\xi\left(t; w_{1}, w_{2}, \alpha\right), \eta\left(t; w_{1}, w_{2}, \alpha\right)\right)$, $(n-1) T \leq t \leq n T$, $n\in \mathbb{N}$ be the order-1 periodic solution. Denote

and define

Theorem 8. For given $w_1$, $w_2$, and $\alpha$ and Case I, $\mathbf{z}(\mathrm{t}) = \left(\xi\left(t; w_{1}, w_{2}, \alpha\right), \eta\left(t; w_{1}, w_{2}, \alpha\right)\right)$, $(n-1) T \leq t \leq n T$, $n\in \mathbb{N}$ is orbitally asymptotically stable if

Proof. For Case I, $\mathbf{z}(\mathrm{t}) = \left(\xi\left(t; w_{1}, w_{2}, \alpha\right), \eta\left(t; w_{1}, w_{2}, \alpha\right)\right)$ lies in the region $S_1$. Then,

so that

Thus, it has

and

Therefore,

To sum up, there is $\mu_{1} < 1$ if

□

Remark 1. For Case Ⅱ, it is only necessary to consider that $\mathbf{\widetilde{z}}(\mathrm{t}) = (\widetilde{\xi}(t), \widetilde{\eta}(t))$ $((n-1) \widetilde{T} \leq t \leq n \widetilde{T})$, $n\in \mathbb{N}$ interacts with $\sum_s$. In such situation, $\mathbf{\widetilde{z}}(\mathrm{t}) = (\widetilde{\xi}(t), \widetilde{\eta}(t))$ $((n-1) \widetilde{T} \leq t \leq n \widetilde{T})$, $n\in \mathbb{N}$ is always orbitally asymptotically stable since it always leaves the sliding region $\sum_s$ from the same point.

4.   Numerical simulations

To illustrate the theoretical results, we present the numerical simulations, which are implemented in MATLAB simulations. The model parameters are assumed to be $K = 150$, $a = 0.3$, $q_1 = 0.02$, $q_2 = 0.02$, $\tilde{r} = (5, 6, 7)$, $\tilde{d} = (0.4, 0.5, 0.6)$, $\tilde{c} = (0.08, 0.1, 0.12)$, which are arbitrarily selected within their reasonable range.

4.1.   Verification of system (2.7)

For system (2.7), three aspects of the imprecision indicator $(w_1, w_2, \alpha)$, the coefficient of refuge $m$, and fishing efforts $\mathbf{E} = (E_1, E_2)$ together affect the dynamical behavior of the system. To illustrate the effect of different parameters on the system, we will verify it by fixing two parameters and changing one parameter.

First, we consider a scenario where fishing activity is not allowed, that is, $\mathbf{E} = \mathbf{0}$. For $m = 0.25$, that is, about 25% prey fishes can go into the shelter, the dynamics of the system (2.7) for different imprecise index on the system are presented in Figures 2–4. It can be observed that for $\alpha = 0.2$, when $w_2$ is small (for example $w_2 = 0.2$), $P^*_1$ is a virtual equilibrium and $P^*_2$ is a real equilibrium; when $w_2$ is big (for example $w_2 = 0.8$), $P^*_1$ is a real equilibrium and $P^*_2$ is a virtual equilibrium; when $w_1$ is small (for example $w_1 = 0.2$), $P^*_1$ is a virtual equilibrium and $P^*_2$ is a real equilibrium for smaller $w_2$; with the increasing of $w_2$ (for example $w_2 = 0.65$), both $P^*_1$ and $P^*_2$ are virtual equilibria; while for bigger $w_2$ (for example $w_2 = 0.8$), $P^*_1$ is a real equilibrium and $P^*_2$ is a virtual equilibrium. Next, for $(w_1, w_2, \alpha) = (0.9, 0.4, 0.9)$, the impact of the refuge coefficient on the dynamics of the system (2.7) is demonstrated, where $m$ is selected to characterize the relative size of the refuge. The dynamics of the system for different $m$ is presented in Figure 5. Obviously, the prey's refuge level does affect the system's dynamic behaviour.

Second, we consider a scenario where fishing activity is allowed. For $m = 0.25$, $(w_1, w_2, \alpha) = (0.5, 0.5, 0.5)$, the dynamics of the system for different $\mathbf{E}$ are presented in Figure 6. It can be observed that the fishing efforts $\mathbf{E}$ have a certain impact on the dynamics of system (2.7). For smaller fishing efforts (for example $\mathbf{E} = (0.1, 0.14)$), $P^*_1$ is a virtual equilibrium and $P^*_2$ is a real equilibrium; with the increasing of fishing efforts (for example $\mathbf{E} = (1, 1.4)$), both $P^*_1$ and $P^*_2$ are virtual equilibria; for bigger fishing efforts (for example $\mathbf{E} = (5, 7)$), $P^*_1$ is a real equilibrium and $P^*_2$ is a virtual equilibrium.

4.2.   Verification of model (2.8)

Given that $m = 0.3$. For control parameters $(w_1, w_2, \alpha) = (0.2, 0.3, 0.5)$, $(E_1, E_2) = (2, 18)$, $l = 1$, when $y_H = 1$, an order-1 periodic solution exists for case I, which totally lies in the region $S_1$, as presented in Figure 7(a); when $y_H = 0.1$, an order-1 periodic solution exists for case Ⅱ, which totally lies in the region $S_2$, as presented in Figure 7(b). While for control parameters $(w_1, w_2, \alpha) = (0.505, 0.505, 0.01)$, $(E_1, E_2) = (2, 18)$, $l = 1$, $y_H = 0.1$, an order-I periodic solution exists for case Ⅱ, which includes a sliding trajectory in region $\sum_s$, as shown in Figure 7(c).

5.   Conclusions

In the natural ecosystem, prey species may exhibit refuge effect when facing the threat of predator species. When the number of predators is relative high compared to the prey, certain percentages of the prey will hide, and when the number of predators is relative small compared to the prey, prey choose not to hide. In addition, natural species are affected by environmental changes in ecosystems, resulting in certain inaccuracies or uncertainties in some key biological parameters. Considering this phenomenon, a Filippov-type fishery model with discontinuous refuge effect and triangle fuzzy number was proposed. Through the Filippov theorem, the sliding mode dynamics of Filippov-type predator-prey system with continuous harvesting were analyzed. The results show that the system may present different types (real, virtual, and pseudo) of equilibrium under different conditions (Theorem 4, Theorem 5, Figures 2–6).

Considering the exploitation of fish resources, a fishing model with threshold control was established by adopting a linear dependent feedback fishing strategy. The complex dynamic properties of the control model are analyzed, including the existence and stability of coexisting order-1 periodic solutions. For two special cases, we provide conditions for the existence of first-order periodic solutions that depend on the relative values of $y_T$ and $n$ (Theorem 6, Theorem 7, Figure 7). The results show that different dynamic behaviors can be obtained in the system with discontinuous refuge and different type of harvesting activities.

Use of AI tools declaration

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

Acknowledgement

The work was supported by the National Natural Science Foundation of China (No. 11401068).

Conflict of interest

The authors declare there is no conflicts of interest.