Review

Cytosine hydroxymethylation by TET enzymes: From the control of gene expression to the regulation of DNA repair mechanisms, and back

  • Received: 06 July 2018 Accepted: 18 September 2018 Published: 25 September 2018
  • Chromatin is a complex multi-scale structure composed of DNA wrapped around nucleosomes. The compaction state is finely regulated mainly by epigenetic marks present not only on nucleosomes but also on the DNA itself. The most studied DNA post-transcriptional modification is 5-methylcytosine (5-mC). Methylation of the cytosine at CpG islands localized at the promoter is associated with repression of transcription. On the contrary, enrichment of 5-hydroxymethylcytosine (5-hmC), one of the oxidation products of 5-mC by TET (ten-eleven translocation) enzymes, on promoters and enhancers promotes transcription activation. Recently, a new role of 5-hmC has been proposed in the context of DNA repair. 5-hmC was found to be enriched at DNA lesions and knockdown of TET led to impaired repair efficiency. Here, we review our current knowledge regarding the role of the regulation of the 5-mC/5-hmC balance by TET enzymes in the context of transcription modulation as well as DNA repair processes. In a final section, we speculate on the potential involvement of TET proteins in DNA repair mechanisms associated with transcription activation.

    Citation: Audrey Lejart, Gilles Salbert, Sébastien Huet. Cytosine hydroxymethylation by TET enzymes: From the control of gene expression to the regulation of DNA repair mechanisms, and back[J]. AIMS Biophysics, 2018, 5(3): 182-193. doi: 10.3934/biophy.2018.3.182

    Related Papers:

    [1] Shuqi Zhai, Qinglong Wang, Ting Yu . Fuzzy optimal harvesting of a prey-predator model in the presence of toxicity with prey refuge under imprecise parameters. Mathematical Biosciences and Engineering, 2022, 19(12): 11983-12012. doi: 10.3934/mbe.2022558
    [2] Parvaiz Ahmad Naik, Muhammad Amer, Rizwan Ahmed, Sania Qureshi, Zhengxin Huang . Stability and bifurcation analysis of a discrete predator-prey system of Ricker type with refuge effect. Mathematical Biosciences and Engineering, 2024, 21(3): 4554-4586. doi: 10.3934/mbe.2024201
    [3] Christian Cortés García . Bifurcations in a discontinuous Leslie-Gower model with harvesting and alternative food for predators and constant prey refuge at low density. Mathematical Biosciences and Engineering, 2022, 19(12): 14029-14055. doi: 10.3934/mbe.2022653
    [4] Yuxuan Zhang, Xinmiao Rong, Jimin Zhang . A diffusive predator-prey system with prey refuge and predator cannibalism. Mathematical Biosciences and Engineering, 2019, 16(3): 1445-1470. doi: 10.3934/mbe.2019070
    [5] Yuanfu Shao . Bifurcations of a delayed predator-prey system with fear, refuge for prey and additional food for predator. Mathematical Biosciences and Engineering, 2023, 20(4): 7429-7452. doi: 10.3934/mbe.2023322
    [6] Rajalakshmi Manoharan, Reenu Rani, Ali Moussaoui . Predator-prey dynamics with refuge, alternate food, and harvesting strategies in a patchy habitat. Mathematical Biosciences and Engineering, 2025, 22(4): 810-845. doi: 10.3934/mbe.2025029
    [7] Jinxing Zhao, Yuanfu Shao . Bifurcations of a prey-predator system with fear, refuge and additional food. Mathematical Biosciences and Engineering, 2023, 20(2): 3700-3720. doi: 10.3934/mbe.2023173
    [8] Xiaoyuan Chang, Junjie Wei . Stability and Hopf bifurcation in a diffusivepredator-prey system incorporating a prey refuge. Mathematical Biosciences and Engineering, 2013, 10(4): 979-996. doi: 10.3934/mbe.2013.10.979
    [9] Christian Cortés García, Jasmidt Vera Cuenca . Impact of alternative food on predator diet in a Leslie-Gower model with prey refuge and Holling Ⅱ functional response. Mathematical Biosciences and Engineering, 2023, 20(8): 13681-13703. doi: 10.3934/mbe.2023610
    [10] Hongqiuxue Wu, Zhong Li, Mengxin He . Dynamic analysis of a Leslie-Gower predator-prey model with the fear effect and nonlinear harvesting. Mathematical Biosciences and Engineering, 2023, 20(10): 18592-18629. doi: 10.3934/mbe.2023825
  • Chromatin is a complex multi-scale structure composed of DNA wrapped around nucleosomes. The compaction state is finely regulated mainly by epigenetic marks present not only on nucleosomes but also on the DNA itself. The most studied DNA post-transcriptional modification is 5-methylcytosine (5-mC). Methylation of the cytosine at CpG islands localized at the promoter is associated with repression of transcription. On the contrary, enrichment of 5-hydroxymethylcytosine (5-hmC), one of the oxidation products of 5-mC by TET (ten-eleven translocation) enzymes, on promoters and enhancers promotes transcription activation. Recently, a new role of 5-hmC has been proposed in the context of DNA repair. 5-hmC was found to be enriched at DNA lesions and knockdown of TET led to impaired repair efficiency. Here, we review our current knowledge regarding the role of the regulation of the 5-mC/5-hmC balance by TET enzymes in the context of transcription modulation as well as DNA repair processes. In a final section, we speculate on the potential involvement of TET proteins in DNA repair mechanisms associated with transcription activation.


    Predator-prey system is one of the most important systems in ecology to study the relationship between prey and predator. For the sake of avoiding predation and increasing the survival rate, prey takes different strategies, for example, changing colour, releasing smell, injecting poison, seeking refuges and so on. Since Gause et al.[1] and Smith [2] introduced a quantity involving prey refuge, the study of prey refuge changing the dynamics of predator-prey system has been widely concerned by researchers, see [3,4,5,6] in details. Olivares and Jiliberto [7] put forward that the quantity on prey refuge can be expressed as a constant number or depending on the density of prey. Kar [8], Huang et al. [9] and Tripathi [10] separately studied the impact of refuge on dynamical behavior of predator-prey systems with Holling type Ⅱ, Holling type Ⅲ and Beddington-DeAngelis type function responses. Han et al. [11] studied the dynamical behavior of a spatiotemporal predator-prey model with refuge and diffusion including boundedness, permanence, coexisting equilibrium and Turing instability parameter space. Qi and Meng [12] investigated the threshold behavior of a stochastic predator-prey system with refuge and fear effect, and obtained the thresholds determining the extinction and persistence in the mean of prey and predator and the existence of a unique ergodic stationary distribution.

    Natural resources are the material basis for human survival and development, among which biological resources are the resources with the most unique development mechanism. They can be continuously renewed via their own reproduction, which is the material embodiment of earth's biological diversity. However, biological resources are also a kind of exhaustible natural resources, and the environmental deterioration or human unreasonable use or excessive development will not only make its quantity and quality decline, but may even lead to the extinction of species. With the demand for food and energy for the development of human beings being rapidly increasing, biological resources are being exploited more and more. Scientific management of the use of marketable biological resources such as fishery, forestry and wildlife is necessary to protect ecosystems. The exploitation of renewable resources used to base on maximum sustainable yield (MSY), but it is impractical to completely ignore the cost operation of resource exploitation. In the face of the insufficiency of MSY, people tried to replace MSY with the optimal sustainable yield (OSY). Lately the research on optimal management of renewable resources has aroused the interest of many scientists and researchers [13,14,15,16,17] and the references therein. Clark [18,19] laid the foundation in the area of work. Kar and Chaudhuri [20] researched a harvesting model on two competing preys and a predator. He and Zhou [21] studied optimal harvesting for a nonlinear hierarchical age-structured population model.

    In the study of biological mathematical models, many works focus on the deterministic model, which is mainly based on the law of large numbers and the assumption that the number of biological individuals is sufficiently large, the behaviors of the system will present a relatively stable statistical regularity, so as to be approximately regarded as a deterministic model. However, any species in nature will inevitably be affected by the complexity of the ecosystem itself and the limitations of human cognition, for example, the uncertainty caused by nature environment or human society such as forest fire, flood disaster, debris flow, volcanic eruption, earthquake, climate warming and so on, the uncertainty caused by simplified hypothesis, the uncertainty caused by replacing infinite samples with limited experimental samples in statistical analysis, the uncertainty caused by test results restricted by objective experimental conditions, even the uncertainty caused by some unknown factors. To tackle these uncertainty factors, many researchers used to adopt interval approach [22,23,24,25,26] and stochastic approach [27,28,29,30,31]. In the interval approach the uncertain parameters are characterized by interval-valued functions. In the stochastic approach the uncertain parameters are replaced as stochastic variables with known probability distribution functions. Nevertheless, the interval approach is too simple to do not achieve ideal effect on complex problems when characterizing imprecise parameters. Also a question arising for stochastic approach is that for imprecise parameters, appropriate probability distribution functions will make the model more complex. To overcome these difficulties, we apply fuzzy set theory to depict imprecise parameters. Fuzzy set theory was first come up with by professor Zadeh [32] and he also considered the appliance of fuzzy differential equations is an inherent way to model dynamic systems under possibilistic uncertainty [33]. And Chang and Zadeh [34] first initiated the concept of the fuzzy derivative. And then Kaleva [35] first put forward the concept of differential equations in a fuzzy environment. All derivatives are considered as Hukuhara derivative or generalized derivative in FDE (fuzzy differential equation). Bede et al. [36] studied first order linear fuzzy differential equations under generalized differentiability, and Khastan and Nieto [37] analysed a boundary value problem for second order fuzzy differential equations. Stefanini and Bede [38] proposed generalization of the Hukuhara difference for compact convex sets, to introduce and study a generalization of the Hukuhara differentiability to the case of interval-valued functions. Puri and Ralescu [39] generalized the concept of Hukuhara differentiability (H-derivative) for set-valued mappings to fuzzy mappings classes. Kaleva [35] developed the theory for FDE by using the H-derivative. Of course, different ways to solve FDE have a huge number of works in this field. Bede and Gal [40] studied the solutions of fuzzy differential equations based on generalized differentiability. Hullermeier [41] proposed an approach to modelling and simulation of uncertain dynamical systems, and showed that all (reasonable) fuzzy functions can be approximated to any degree of accuracy in this way. Allahviranloo and Ahmadi [42] applied Laplace transform method to solve linear FDE. Also the existence and uniqueness of solution for fuzzy initial value problems are investigated by Villamizar-Roa et al. [43]. Different forms of FDEs are come up with in dynamics and good results have been achieved [23,44,45,46,47]. Bassanezi et al. [48] is the cornerstone of studying the stability of dynamical system by employing fuzzy differential equations. Mizukoshi et al. [49] studied the stability of fuzzy dynamical systems considering initial conditions under fuzziness. Guo et al. [50] established fuzzy impulsive functional differential equations and applied them in logistic model and Gompertz model. The existence and uniqueness of solution for fuzzy functional differential equations were obtained by Lupulescu [51]. Recently, Sadhukhan et al. [52] considered optimal harvesting of a food chain model in fuzzy environment, it is worth mentioning that based on the fuzzy instantaneous annual rate of discount, the optimal harvesting policy is discussed. Pal and Mahapatra [23] drew interval number and fuzzy number into prey and predator interaction's model. The existence and stability of biological and bionomic equilibria are investigated under interval parameters, and the optimal harvesting policy considering instantaneous annual rate of discount in fuzzy environment is studied. Besides, Pal et al. [44] discussed the stability of a fuzzy predator-prey harvesting model with toxicity. And Pal et al. [53] also put forward a predator-prey model with fuzzy optimal harvesting in the environment of toxicity. The former is a major consideration on fuzzy number applied in the parameters of model, the latter mainly considers the inflation and discount rates as fuzzy number, and obtains the fuzzy optimal harvesting. However, for all we know, no one tried to consider both fuzzy biological parameters and fuzzy optimal harvesting in a model. This paper first establishes a predator-prey model with refuge, simultaneously considers both fuzzy biological parameters and fuzzy optimal harvesting into model. The dynamical behaviors covering the existence and stability of biological equilibria are performed. Moreover, there exist bionomic equilibria under fuzzy biological parameters. The optimal harvesting policy is obtained under imprecise inflation and discount in fuzzy environment.

    The remaining part of this article is emerged as follows. Firstly, the concepts on fuzzy set and weighted sum method are put forward in Section 2. Section 3 presents the suitable model under fuzzy biological parameters in imprecise environment. The positivity and boundedness of system, as well as the existence and stability of equilibria are investigated in Sections 4 and 5, respectively. In Section 6, we discuss all kinds of bionomic equilibria. The fuzzy optimal harvesting is presented in Section 7. Some numerical examples verifying the theoretical results are displayed in Section 8. The last section gives a brief summary of the work of this paper.

    This section provides definitions, lemmas and methods that will be used in latter sections.

    Definition 2.1. [32] Fuzzy set: A fuzzy set $ \tilde{A} $ in a universe of discourse $ X $ is defined as the following set of pairs $ \tilde{A} = \{(x, \mu_{\tilde{A}}(x)):x\in X\} $. The mapping $ \mu_{\tilde{A}}:X\rightarrow [0, 1] $ is called the membership function of the fuzzy set $ \tilde{A} $ and $ \mu_{\tilde{A}} $ is called the membership value or degree of membership of $ x\in X $ in the fuzzy set $ \tilde{A} $.

    Definition 2.2. [54] $ \alpha $-cut of fuzzy set: The $ \alpha $-cut of a fuzzy set $ \tilde{A} $ is a crisp set and it is defined by $ A_\alpha = \{x:\mu_{\tilde{A}}(x)\geq\alpha \}, \alpha\in(0, 1] $. For $ \alpha = 0 $ the support of $ \tilde{A} $ is defined as $ A_0 = Supp(\tilde{A}) = \overline{\{x\in R, \mu_{\tilde{A}}(x) > 0\}} $.

    Definition 2.3. [32] Convex fuzzy set: A convex fuzzy set $ \tilde{A} $ is a fuzzy set on a continuous universe such that for all $ \alpha $, $ A_\alpha $ is a convex classical set.

    Definition 2.4. [55] Fuzzy number: A fuzzy number is a convex fuzzy set with $ X = R $.

    Definition 2.5. [56] Triangular fuzzy number: A triangular fuzzy number (TFN) $ \tilde{A}\equiv(a_1, a_2, a_3) $ is fuzzy set of the real line $ R $ characterized by the membership function $ \mu_{\tilde{A}}:R\rightarrow [0, 1] $ as follows

    $ \mu_{\tilde{A}} = \left\{ xa1a2a1  if  a1xa2,a3xa3a2  if  a2xa3,0  otherwise.
    \right. $

    Therefore, the $ \alpha $-cut of triangular fuzzy number is a bounded and closed internal $ [A_L(\alpha), A_R(\alpha)] $, where $ A_L(\alpha) = {\inf}\{x:\mu_{\tilde{A}}(x)\geq\alpha\} = a_1+\alpha(a_2-a_1) $ and $ A_R(\alpha) = {\sup}\{x:\mu_{\tilde{A}}(x)\geq\alpha\} = a_3-\alpha(a_3-a_2) $.

    The arithmetic operations on fuzzy numbers are provided by the following lemmas. Firstly, let us define $ E $ as the set of all upper semi-continuous normal convex fuzzy numbers with $ \tilde{A}\in E $, then the $ \alpha $-level set $ A_\alpha $ is a bounded closed interval which can be expressed as $ A_\alpha = [A_L^\alpha, A_R^\alpha] $, where $ A_L^\alpha = \inf\{x:\mu_{\tilde{A}}(x)\geq\alpha\} $ and $ A_R^\alpha = \sup\{x:\mu_{\tilde{A}}(x)\geq\alpha\} $.

    Lemma 2.1. [57] If $ \tilde{A}, \tilde{B}\in E $, then for $ \alpha\in(0, 1] $,

    $ \rm(1) $ $ [\tilde{A}+\tilde{B}]_\alpha = [A_L^\alpha+B_L^\alpha, A_R^\alpha+B_R^\alpha] $,

    $ \rm(2) $ $ [\tilde{A}\cdot\tilde{B}]_\alpha = [\min\{A_L^\alpha B_L^\alpha, A_L^\alpha B_R^\alpha, A_R^\alpha B_L^\alpha, A_R^\alpha B_R^\alpha\}, \max\{A_L^\alpha B_L^\alpha, A_L^\alpha B_R^\alpha, A_R^\alpha B_L^\alpha, A_R^\alpha B_R^\alpha\}] $,

    $ \rm(3) $ $ [\tilde{A}-\tilde{B}]_\alpha = [A_L^\alpha-B_R^\alpha, A_R^\alpha-B_L^\alpha] $.

    Lemma 2.2. [58] $ [A_L^\alpha, A_R^\alpha], 0 < \alpha\leq1 $ is a family of nonempty intervals. If

    $ (1) $ $ [A_L^\alpha, A_R^\alpha]\supset[A_L^\beta, A_R^\beta] $ for $ 0 < \alpha\leq\beta $ and

    $ (2) $ $ [\lim_{k\rightarrow \infty}A_L^{\alpha_k}, \lim_{k\rightarrow \infty}A_R^{\alpha_k}] = [A_L^\alpha, A_R^\alpha] $,

    whenever non-decreasing sequence $ \alpha_k $ converges to $ \alpha\in(0, 1] $; then, the family $ [A_L^\alpha, A_R^\alpha] $, $ 0 < \alpha\leq1 $, represents the $ \alpha $-level sets of a fuzzy number $ \tilde{A} $ in $ E $. Conversely, if $ [A_L^\alpha, A_R^\alpha] $, $ 0 < \alpha\leq1 $ are the $ \alpha $-level sets of a fuzzy number $ \tilde{A}\in E $, then the conditions ${\rm(1)}$ and ${\rm(2)}$ hold true.

    In weighted sum method [59], a simple utility function is in the form of $ w_ig_i $ for $ i $th objective, where $ w_i $ stands for the weight of $ i $th objective. And the total utility function $ U $ is expressed as follows

    $ U=liwigi,   i=1,2,,l,
    $
    (2.1)

    where $ w_i > 0 $ and $ \sum_i^l w_i = 1 $ are satisfied.

    Motivated by the predator-prey system (3) in Wang et al. [5] as follows

    $ {dxdt=rx(1xK)c(1m)xyq1E1x,dydt=dysy2+e(1m)xyq2E2y,
    $
    (3.1)

    where $ r, K, c, m, d, s $ and $ e $ are positive constants. $ x(t) $ and $ y(t) $ depict the scale of prey and predator at $ t $ with initial densities $ x_1(0) > 0 $ and $ x_2(0) > 0 $. $ r $ is the growth rate of prey, $ K $ stands for the carrying capacity of prey, $ c $ is the predation rate for prey, $ m $ shows the proportion of quantity of prey using refuges, $ d $ expresses as the death rate of predator and $ s $ is the intraspecies competition rate of predator, $ c $ is the predation rate, $ e $ represents the product of predation rate and conversion rate, $ q_1, q_2 $ and $ E_1, E_2 $ respectively say the catchability coefficients and the harvesting efforts of prey and predator.

    On account of the concept of fuzzy set, we regard the imprecise parameters $ \tilde{r}, \tilde{c}, \tilde{d}, \tilde{s} $ and $ \tilde{e} $ as triangular fuzzy numbers. Throughout Hukuhara's concept of derivative [39], we present fuzzy differential equations to describe predator-prey interaction with harvesting items as follows

    $ {~dxdt=˜rx(1xK)˜c(1m)xyq1E1x,~dydt=˜dy˜sy2+˜e(1m)xyq2E2y.
    $
    (3.2)

    Here, prey increases in the absence of predator with fuzzy growth rate $ \tilde{r} $ while predator decreases with fuzzy mortality $ \tilde{d} $ for lack of prey. $ \tilde{s} $ expresses as fuzzy decrease rate of predator attributed to intraspecific competition. The fuzzy quantity of the prey predated per one predator under unit time yields to $ \tilde{c}(1-m)x $, and $ \tilde{e} $ represents the fuzzy growth rate of predator attacking on prey.

    Applying $ \alpha $-level to cut triangular fuzzy numbers $ \tilde{r}, \tilde{c}, \tilde{d}, \tilde{s} $ and $ \tilde{e} $, system (3.2) can be expressed as follows

    $ {(dxdt)αL=(rL)αx(rR)αx2K(cR)α(1m)xyq1E1x,(dxdt)αR=(rR)αx(rL)αx2K(cL)α(1m)xyq1E1x,(dydt)αL=(dR)αy(sR)αy2+(eL)α(1m)xyq2E2y,(dydt)αR=(dL)αy(sL)αy2+(eR)α(1m)xyq2E2y.
    $
    (3.3)

    By the method of weighted sum, the above system translates into

    $ {dxdt=w1(dxdt)αL+w2(dxdt)αR,dydt=w1(dydt)αL+w2(dydt)αR,
    $
    (3.4)

    where $ w_1 $ and $ w_2 $ express as weight satisfying $ w_1+w_2 = 1 $, and $ w_1, w_2\geq0 $. The simplified form of system (3.4) is given by

    $ {dxdt=A1xA2x2KA3(1m)xyq1E1x,dydt=B1yB2y2+B3(1m)xyq2E2y,
    $
    (3.5)

    where

    $ A1=w1(rL)α+w2(rR)α,  A2=w1(rR)α+w2(rL)α,  A3=w1(cR)α+w2(cL)α,B1=w1(dR)α+w2(dL)α,  B2=w1(sR)α+w2(sL)α,  B3=w1(eL)α+w2(eR)α.
    $
    (3.6)

    Of course if we consider imprecise biological parameters with interval parameters rather than fuzzy parameters, then system (3.5) can be translated into system (5) in Wang [5]. Moreover, if we neglect the intraspecific competition rates of prey and predator, as well as refuges, system (3.5) can be simplified to Eqs (7) and (8) in Pal et al. [22].

    Throughout the following theorem, in this section, we make sure the positivity and boundedness of system (3.5).

    Theorem 4.1. Any solution of system ${\rm(3.5)}$ is positive and bounded satisfying initial conditions $ x(0) > 0 $ and $ y(0) > 0 $ for all $ t > 0 $.

    Proof. The right hand of system $(3.5)$ satisfies completely continuous and locally Lipschitzian on $ C $, which guarantees that there exists a unique solution $ (x(t), y(t)) $ of system $(3.5)$ under initial conditions $ x(0) > 0 $ and $ y(0) > 0 $ on $ [0, \xi) $, where $ 0 < \xi < +\infty $. It follows from system $(3.5)$ with $ x(0) > 0 $ and $ y(0) > 0 $ that

    $ x(t)=x(0)exp[t0(A1A2x(s)KA3(1m)y(s)q1E1)ds]>0,y(t)=y(0)exp[t0(B1B2y(s)+B3(1m)x(s)q2E2)ds]>0,
    $
    (4.1)

    which ensures the solution of system (3.5) is positive.

    On the other side, the first differential equation of system (3.5) provides

    $ dxdtA1x(1xA1KA2),
    $
    (4.2)

    which implies $ \limsup_{t\rightarrow \infty}x(t)\leq\frac{A_1K}{A_2}\equiv k_1 $. From the second differential equation of system (3.5), one has

    $ \frac{dy}{dt}\leq k_1B_3(1-m)y-B_2y^2, $

    which can be written as

    $ dydtk1B3(1m)y(1yk1B3(1m)B2).
    $
    (4.3)

    It is easy to see $ \limsup_{t\rightarrow \infty}y(t)\leq\frac{k_1B_3(1-m)}{B_2} $. Therefore, the solution of system (3.5) is bounded.

    We firstly put forward the existence of equilibria of system (3.5). Besides the local stability and global stability of the above equilibria are investigated by the method of Jacobian matrix and Lyapunov function, respectively.

    By a simple calculation, the equilibria of system (3.5) are given by

    (1) Trivial equilibrium: $ S_0 = (0, 0) $.

    (2) Axial equilibrium: $ S_1 = (\frac{K(A_1-q_1E_1)}{A_2}, 0) $ exists if $ A_1 > q_1E_1 $.

    (3) Axial equilibrium: $ S_2 = (0, -\frac{q_2E_2+B_1}{B_2}) $ does not exist since $ -\frac{q_2E_2+B_1}{B_2} < 0 $.

    (4) Interior equilibrium: $ S^* = (x^*, y^*) $, where

    $ x=KA1B2+KA3B1(1m)+KA3(1m)q2E2q1E1B2KA2B2+KA3B3(1m)2,y=B3(1m)x(B1+q2E2)B2
    $
    (5.1)

    exists if $ B_3(1-m)x^* > B_1+q_2E_2 $ and $ A_1B_2+A_3(1-m)(B_1+q_2E_2) > B_2q_1E_1 $ are satisfied.

    On account of the following theorem, the local stability of equilibria $ S_0, S_1, S^* $ for system (3.5) is discussed.

    Theorem 5.1. The following conclusions hold:

    $ \rm(1) $ Trivial equilibrium $ S_0 $ becomes stable node when $ A_1 < q_1E_1 $ or saddle point when $ A_1 > q_1E_1 $.

    $ \rm(2) $ Axial equilibrium $ S_1 $ becomes stable node when

    $ A1>q1E1  and  A2(B1+q2E2)>KB3(1m)(A1q1E1),
    $
    (5.2)

    or saddle point when

    $ A1>q1E1  and  A2(B1+q2E2)<KB3(1m)(A1q1E1).
    $
    (5.3)

    $ \rm(3) $ Interior equilibrium $ S^* $ is stable node or focus when

    $ B3(1m)x>B1+q2E2  and  A1B2+A3(1m)(B1+q2E2)>B2q1E1.
    $
    (5.4)

    Proof. The following Jacobian matrix of system (3.5) is expressed as

    $ M=(A1q1E12A2xKA3(1m)yA3(1m)xB3(1m)y(B1+q2E2)2B2y+B3(1m)x).
    $
    (5.5)

    (1) The Jacobian matrix $ M_0 $ at $ S_0 $ is represented as

    $ M_0 = (A1q1E100B1q2E2)
    . $

    According to the characteristic equation about the above matrix

    $ (\lambda-A_1+q_1E_1)(\lambda+B_1+q_2E_2) = 0, $

    one has $ \lambda_1 = A_1-q_1E_1, \ \ \lambda_2 = -B_1-q_2E_2 $, hence $ A_0 $ is stable node when $ A_1 < q_1E_1 $ or saddle point when $ A_1 > q_1E_1 $.

    (2) The Jacobian matrix $ M_1 $ at $ S_1 $ is depicted as

    $ M_1 = ((A1q1E1)KA3(1m)(A1q1E1)A20(B1+q2E2)+KB3(1m)(A1q1E1)A2)
    . $

    In accordance with the characteristic equation of the above matrix

    $ (\lambda+A_1-q_1E_1)\Bigl(\lambda+B_1+q_2E_2-\frac{KB_3(1-m)(A_1-q_1E_1)}{A_2}\Bigr) = 0, $

    we have $ \lambda_1 = -(A_1-q_1E_1) $ and $ \lambda_2 = -(B_1+q_2E_2)+\frac{KB_3(1-m)(A_1-q_1E_1)}{A_2} $. Obviously, $ S_1 $ exists if $ A_1 > q_1E_1 $ holds, hence $ S_1 $ is a stable node when $ A_1 > q_1E_1 $ and $ A_2(B_1+q_2E_2) > KB_3(1-m)(A_1-q_1E_1) $ or saddle point when $ A_1 > q_1E_1 $ and $ A_2(B_1+q_2E_2) < KB_3(1-m)(A_1-q_1E_1) $.

    (3) The Jacobian matrix $ M^* $ at $ S^* $ is written as

    $ M^* = (A2xKA3(1m)xB3(1m)yB2y)
    . $

    The characteristic equation is expressed as

    $ |λ+A2xKA3(1m)xB3(1m)yλ+B2y|=λ2+(B2y+A2xK)λ+A3B3(1m)2xy=0.
    $
    (5.6)

    It is easy to see the sum of the roots of Eq (5.6) is negative, i.e., $ \lambda_1+\lambda_2 = -\Bigl(B_2y^*+\frac{A_2x^*}{K}\Bigr) < 0 $ and the product of the roots of Eq (5.6) is positive, i.e., $ \lambda_1\lambda_2 = A_3B_3(1-m)^2x^*y^* > 0 $. Thus the roots of quadratic Eq (5.6) are negative real number or complex conjugates with negative real parts. Therefore, $ (x^*, y^*) $ is stable node or focus when the existence conditions $ B_3(1-m)x^* > B_1+q_2E_2 $ and $ A_1B_2+A_3(1-m)(B_1+q_2E_2) > B_2q_1E_1 $ are satisfied.

    We next discuss the global stability of interior equilibrium $ S^* $.

    Theorem 5.2. If interior equilibrium $ S^* $ exists, then it is globally asymptotically stable.

    Proof. The Lyapunov function is established as follows

    $ V(x,y)=xxxln(xx)+l[yyyln(yy)],
    $
    (5.7)

    where $ l > 0 $ is a constant. Differentiating both sides of the above equation in regard to time $ t $, one has

    $ dVdt=xxxdxdt+lyyydydt=(xx)[A1A2xKA3(1m)yq1E1]+l(yy)[B1B2y+B3(1m)xq2E2]=A2K(xx)2A3(1m)(xx)(yy)lB2(yy)2+lB3(1m)(xx)(yy).
    $

    Choose $ l = \frac{A_3}{B_3} $, it yields to

    $ \frac{dV}{dt} = -\frac{A_2}{K}(x-x^*)^2-\frac{A_3B_2}{B_3}(y-y^*)^2 < 0. $

    Since $ \frac{dV}{dt} $ in some neighborhood of $ (x^*, y^*) $ is negative semidefinite, the equilibrium $ (x^*, y^*) $ is globally asymptotically stable.

    As we have discussed in Section 5, the biological equilibrium yields to $ \frac{dx}{dt}\equiv0 $ and $ \frac{dy}{dt}\equiv0 $. And economic equilibrium is expressed by the total revenue equaling the total cost. Then bionomic equilibrium is a combination of biological equilibrium and economic equilibrium. In this section, we analyse the existence of all bionomic equilibria.

    Firstly, $ c_1, c_2 $ denote unit capture cost for prey $ x $ and predator $ y $, respectively. And $ p_1, p_2 $ show unit biomass price for prey $ x $ and predator $ y $, respectively. Define the profit function as follows

    $ N=(p1q1xc1)E1+(p2q2yc2)E2=N1+N2,
    $
    (6.1)

    where $ N_1 = (p_1q_1x-c_1)E_1, N_2 = (p_2q_2y-c_2)E_2 $. $ N_1 $ and $ N_2 $ separately represent the profit function of prey $ x $ and predator $ y $. And the bionomic equilibrium $ (x_\infty, y_\infty, E_{1\infty}, E_{2\infty}) $ is governed by

    $ {A1xA2x2KA3(1m)xyq1E1x=0,B1yB2y2+B3(1m)xyq2E2y=0,(p1q1xc1)E1+(p2q2yc2)E2=0.
    $
    (6.2)

    The following we consider four cases to excavate the bionomic equilibria of system (3.5).

    Case Ⅰ Boundary bionomic equilibria with no catching species $ x $: $ E_{1\infty} = 0 $. One has

    $ {p2q2yc2=0,A1xA2x2KA3(1m)xy=0,B1B2y+B3(1m)xq2E2=0.
    $
    (6.3)

    Clearly, $ (x_\infty, y_\infty, E_{1\infty}, E_{2\infty}) = (0, \frac{c_2}{p_2q_2}, 0, -\frac{p_2q_2B_1+B_2c_2}{p_2q_2^2}) $ satisfies Eq (6.3), we have to omit the bionomic equilibrium owing to $ E_{2\infty} < 0 $.

    Then solving the second equation for a nonzero $ x_\infty $ yields

    $ x=A1Kp2q2A3(1m)c2KA2p2q2,
    $
    (6.4)

    which exists if satisfying $ A_1p_2q_2 > A_3(1-m)c_2 $. Substituting Eq (6.4) into the third equation of Eq (6.3) gives

    $ E2=A1B3Kp2q2(1m)A3B3c2K(1m)2A2B1p2q2A2B2c2A2p2q22,
    $
    (6.5)

    which is positive provided that

    $ A1B3Kp2q2(1m)>A3B3c2K(1m)2+A2B1p2q2+A2B2c2
    $
    (6.6)

    holds. Therefore, there exists a boundary bionomic equilibrium $ (x_\infty, y_\infty, 0, E_{2\infty}) $.

    Case Ⅱ Boundary bionomic equilibria with no catching species $ y $: $ E_{2\infty} = 0 $. It follows from

    $ {p1q1xc1=0,A1A2xKA3(1m)yq1E1=0,B1yB2y2+B3(1m)xy=0,
    $
    (6.7)

    that $ (x_\infty, y_\infty, E_{1\infty}, E_{2\infty}) = (\frac{c_1}{p_1q_1}, 0, \frac{A_1p_1q_1K-A_2c_1}{p_1q_1^2K}, 0) $ and we get a nonegative semi-trivial bionomic equilibrium $ (x_\infty, 0, E_{1\infty}, 0) $ with condition $ A_1p_1q_1K > A_2c_1 $.

    Then solving the third equation for a nonzero $ y_\infty $ yields

    $ y=B3c1(1m)B1p1q1B2p1q1>0
    $
    (6.8)

    satisfying $ B_3c_1(1-m) > B_1p_1q_1 $. Substituting Eq (6.8) into the second equation of Eq (6.7) gives

    $ E1=KA1B2p1q1+KA3B1(1m)p1q1A2B2c1KA3B3c1(1m)2KB2p1q21,
    $
    (6.9)

    which is positive provided that

    $ KA1B2p1q1+KA3B1(1m)p1q1>A2B2c1+KA3B3c1(1m)2
    $
    (6.10)

    holds. Thus, there exists a boundary bionomic equilibrium $ (x_\infty, y_\infty, E_{1\infty}, 0) $.

    Case Ⅲ Boundary bionomic equilibria with no catching species $ x $ and $ y $: $ E_{1\infty} = 0 $ and $ E_{2\infty} = 0 $.

    In this case, the cost exceeds revenue for both species $ x $ and $ y $, we have to stop harvesting both of them. The bionomic equilibrium is equal to the biological equilibrium with $ E_{1\infty} = E_{2\infty} = 0 $.

    Case Ⅳ Positive bionomic equilibrium with effort fishing on both species $ x $ and $ y $: $ E_{1\infty} > 0 $ and $ E_{2\infty} > 0 $.

    Solving the third equation of Eq (6.2), we obtain $ x_\infty = \frac{c_1}{p_1q_1}, y_\infty = \frac{c_2}{p_2q_2} $. Then, it follows from the second and third equations of Eq (6.2) that

    $ E1=KA1p1p2q1q2A2c1p2q2KA3c2(1m)p1q1Kp1p2q21q2,E2=B3c1(1m)p2q2B1p1p2q1q2B2c2p1q1p1p2q1q22,
    $
    (6.11)

    which exist if satisfying the following conditions

    $ {KA1p1p2q1q2>A2c1p2q2+KA3c2(1m)p1q1,B3c1(1m)p2q2>B1p1p2q1q2+B2c2p1q1.
    $
    (6.12)

    Therefore, there exists the positive bionomic equilibrium $ (x_\infty, y_\infty, E_{1\infty}, E_{2\infty}) $.

    Summarizing the above, we put forward the theorem as follows.

    Theorem 6.1. The following conclusions hold:

    $ \rm(1) $ The boundary bionomic equilibrium $ (x_\infty, y_\infty, 0, E_{2\infty}) $ exists if $ A_1p_2q_2 > A_3(1-m)c_2 $ and Eq $(6.6) $ are satisfied, where $ y_\infty = \frac{c_2}{p_2q_2} $ and $ x_\infty $ and $ E_{2\infty} $ are defined in ${Eqs}$ $ \rm(6.4) $ and ${\rm(6.5)}$, respectively.

    $ \rm(2) $ The boundary bionomic equilibrium $ (x_\infty, 0, E_{1\infty}, 0) $ exists if $ A_1p_1q_1K>A_2c_1 $ holds, where $ x_\infty = \frac{c_1}{p_1q_1} $ and $ E_{1\infty} = \frac{A_1p_1q_1K-A_2c_1}{p_1q_1^2K} $.

    $ \rm(3) $ The boundary bionomic equilibrium $ (x_\infty, y_\infty, E_{1\infty}, 0) $ exists if $ B_3c_1(1-m)>B_1p_1q_1 $ and ${Eq}$ $(6.10)$ are satisfied, where $ x_\infty = \frac{c_1}{p_1q_1} $ and $ y_\infty $ and $ E_{1\infty} $ are defined in ${Eqs}$ ${\rm(6.8)}$ and ${\rm(6.9)}$, respectively.

    $ \rm(4) $ The boundary bionomic equilibrium $ (x_\infty, y_\infty, 0, 0) $ is equal to the biological equilibrium with $ E_{1\infty} = E_{2\infty} = 0 $ described in Section ${\rm5}$.

    $ \rm(5) $ The positive bionomic equilibrium $ (x_\infty, y_\infty, E_{1\infty}, E_{2\infty}) $ exists if ${Eq}$ ${\rm(6.12)}$ holds, where $ x_\infty = \frac{c_1}{p_1q_1} $, $ y_\infty = \frac{c_2}{p_2q_2} $ and $ E_{1\infty} $ and $ E_{2\infty} $ are defined in ${Eq}$ ${\rm(6.11)}$.

    This section studies the optimal harvesting strategy with fuzzy net discount rate of inflation. Previously, denote $ \tilde{k} $ and $ \tilde{r} $ the inflation and discount rates and they are considered as fuzzy number essentially. The net discount rate of inflation $ \tilde{\delta} $ is the difference of the above two fuzzy numbers, that is $ \tilde{\delta} = \tilde{r}-\tilde{k} $, and it can be regarded as triangular fuzzy number, that is $ \tilde{\delta} = (\delta_1, \delta_2, \delta_3) $. Our target is to maximize the value $ \tilde{J} $ as follows

    $ ˜J(E1,E2)=0e˜δt[(p1q1xc1)E1(t)+(p2q2yc2)E2(t)]dt.
    $
    (7.1)

    Using Pontryagin's maximal principle [60] to solve the above fuzzy optimization problem, the harvesting strategies not only guarantee profit maximization, but also maintain an optimal level for species. The control variable $ E_i(t) $ subject to $ 0\leq E_i(t)\leq E_i^{\max}, i = 1, 2 $. According to Maity and Maiti [61], Sadhukhan et al. [52] and Pal and Mahapatra [23], using $ \alpha' $ to cut triangular fuzzy number $ \tilde{\delta} = (\delta_1, \delta_2, \delta_3) $, that is it regarded as an interval number $ [\delta_L, \delta_R] $, the optimization problem Maximize $ \tilde{J}(E_1, E_2) $ is equal to Maximize $ [J_L(E_1, E_2), J_R(E_1, E_2)] $, where

    $ JL(E1,E2)=0eδRt[(p1q1xc1)E1(t)+(p2q2yc2)E2(t)]dt,JR(E1,E2)=0eδLt[(p1q1xc1)E1(t)+(p2q2yc2)E2(t)]dt,δL=δ1+α(δ2δ1),   δR=δ3α(δ3δ2).
    $
    (7.2)

    Again, by the method of weighted sum, it yields to

    $ Max˜J=Max[JL,JR]=Max(w1JL+w2JR),
    $
    (7.3)

    where $ w_1', w_2'\geq0 $ are two weights such that $ w_1'+w_2' = 1 $. The Hamiltonian can be written as

    $ H=(w1eδRt+w2eδLt)[(p1q1xc1)E1+(p2q2yc2)E2] +λ1[A1xA2x2KA3(1m)xyq1E1x] +λ2[B1yB2y2+B3(1m)xyq2E2y],
    $
    (7.4)

    where $ \lambda_i = \lambda_i(t)(i = 1, 2) $ are the adjoint variables. Applying Pontryagin's maximum principle [60], we get the adjoint equations

    $ dλ1dt=Hx,dλ2dt=Hy.
    $
    (7.5)

    It follows from Eqs (7.4) and (7.5) that

    $ dλ1dt={p1q1(w1eδRt+w2eδLt)E1+λ2B3(1m)y   +λ1[A12A2KxA3(1m)yq1E1]},dλ2dt={p2q2(w1eδRt+w2eδLt)E2λ1A3(1m)x   +λ2[B12B2y+B3(1m)xq2E2]}.
    $
    (7.6)

    Substituting the interior equilibrium into Eq (7.6) yields

    $ dλ1dt=λ1A2Kxλ2B3(1m)yp1q1(w1eδRt+w2eδLt)E1,dλ2dt=λ1A3(1m)x+λ2B2yp2q2(w1eδRt+w2eδLt)E2.
    $
    (7.7)

    Now, we neglect $ \lambda_1 $ in Eq (7.7) to get a reduced differential equation for $ \lambda_2 $

    $ a0d2λ2dt2+a1dλ2dt+a2λ2=M2LeδRt+M2ReδLt,
    $
    (7.8)

    where

    $ a0=1,  a1=(A2Kx+B2y),  a2=[A3B3(1m)2+A2B2K]xy,M2L=A3(1m)p1q1xw1E1+A2Kp2q2xw1E2+p2q2w1δRE2,M2R=A3(1m)p1q1xw2E1+A2Kp2q2xw2E2+p2q2w2δLE2.
    $

    The solution of Eq (7.8) as follows

    $ λ2=C1em1t+C2em2t+M2LNLeδRt+M2RNReδLt,
    $
    (7.9)

    where $ C_i(i = 1, 2) $ are arbitrary constants and $ m_i(i = 1, 2) $ are the roots of the auxiliary equations $ a_0m^2+a_1m+a_2 = 0 $ and

    $ N_L = a_0\delta_R^2-a_1\delta_R+a_2\neq0, \ \ \ N_R = a_0\delta_L^2-a_1\delta_L+a_2\neq0. $

    From Eq (7.9), $ \lambda_2 $ is bounded if and only if $ m_i < 0 (i = 1, 2) $ or $ C_i = 0 $. Since it is difficult to check whether $ m_i < 0 $, we consider $ C_i = 0 $. Eq (7.9) is simplified as

    $ λ2=M2LNLeδRt+M2RNReδLt.
    $
    (7.10)

    Similarly, we get

    $ λ1=M1LNLeδRt+M1RNReδLt,
    $
    (7.11)

    where

    $ M1L=B3(1m)yp2q2w1E2+p1q1w1δRE1+B2yp1q1w1E1,M1R=B3(1m)yp2q2w2E2+p1q1w2δLE1+B2yp1q1w2E1.
    $
    (7.12)

    It is easy to see the shadow prices $ \lambda_ie^{\delta_Lt}(i = 1, 2) $ of the two species satisfy the transversality condition at $ \infty $ [18], i.e., they are bounded as $ t\rightarrow \infty $. The Hamiltonian in Eq (7.4) must be maximized for $ E_i\in[0, E_i^{\max}] $. Suppose that the optimal equilibrium does not appear at either $ E_i = 0 $ or $ E_i = E_i^{\max} $, we therefore consider the singular controls

    $ HE1=(w1eδRt+w2eδLt)(p1q1xc1)λ1q1x=0,HE2=(w1eδRt+w2eδLt)(p2q2yc2)λ2q2y=0,
    $
    (7.13)

    that is,

    $ (w1eδRt+w2eδLt)(p1q1xc1)=λ1q1x,(w1eδRt+w2eδLt)(p2q2yc2)=λ2q2y.
    $
    (7.14)

    We substitute $ \lambda_1 $ and $ \lambda_2 $ in Eq (7.11) and Eq (7.10) into Eq (7.14) and obtain

    $ PLeδRt+PReδLt=(w1eδRt+w2eδLt)c1,QLeδRt+QReδLt=(w1eδRt+w2eδLt)c2,
    $
    (7.15)

    where

    $ PL=(p1w1M1LNL)q1x,   PR=(p1w2M1RNR)q1x,QL=(p2w1M2LNL)q2y,   QR=(p2w2M2RNR)q2y.
    $
    (7.16)

    Equation (7.15) is equivalent to the following equations

    $ (PLc1w1)e(δRδL)t=(PRc1w2),(QLc2w1)e(δRδL)t=(QRc2w2).
    $
    (7.17)

    To eliminate $ e^{-(\delta_R-\delta_L)t} $ in Eq (7.17), we consider the following two cases:

    Case Ⅰ If $ P_L\neq c_1w_1' $ and $ Q_L\neq c_2w_1' $, then Eq (7.17) is equal to

    $ (PLc1w1)(QRc2w2)=(QLc2w1)(PRc1w2).
    $
    (7.18)

    Case Ⅱ If $ P_L = c_1w_1' $ or $ Q_L = c_2w_1' $, we have

    $ P_L-c_1w_1' = P_R-c_1w_2' = 0, \ \ \mbox{or}\ \ Q_L-c_2w_1' = Q_R-c_2w_2' = 0. $ (7.18')

    Motivated by the method in [25], firstly transforming Eq (7.14) into the following form

    $ λ1=(w1eδRt+w2eδLt)(p1c1q1x),λ2=(w1eδRt+w2eδLt)(p2c2q2y).
    $
    (7.19)

    Differentiating $ \lambda_i(i = 1, 2) $ in Eq (7.19) with respect to $ t $, it derives that

    $ dλ1dt=(δRw1eδRtδLw2eδLt)(p1c1q1x)+c1q1x(w1eδRt+w2eδLt)[A1A2KxA3(1m)yq1E1],dλ2dt=(δRw1eδRtδLw2eδLt)(p2c2q2y)+c2q2y(w1eδRt+w2eδLt)[B1B2y+B3(1m)xq2E2].
    $
    (7.20)

    Substituting $ \lambda_i(i = 1, 2) $ in Eq (7.19) into Eq (7.6) yields

    $ dλ1dt=p1q1E1(w1eδRt+w2eδLt)B3(1m)y(p2c2q2y)(w1eδRt+w2eδLt)(p1c1q1x)[A12A2KxA3(1m)yq1E1](w1eδRt+w2eδLt),dλ2dt=p2q2E2(w1eδRt+w2eδLt)+A3(1m)x(p1c1q1x)(w1eδRt+w2eδLt)(p2c2q2y)[B12B2y+B3(1m)xq2E2](w1eδRt+w2eδLt).
    $
    (7.21)

    Combined Eqs (7.20) and (7.21) with $ E_1 $ and $ E_2 $ omitted, one has

    $ (p1c1q1x)(δRw1eδRt+δLw2eδLt)=B3(1m)y(p2c2q2y)(w1eδRt+w2eδLt)  +p1[A12A2KxA3(1m)y](w1eδRt+w2eδLt)  +c1A2q1K(w1eδRt+w2eδLt),(p2c2q2y)(δRw1eδRt+δLw2eδLt)=A3(1m)x(p1c1q1x)(w1eδRt+w2eδLt)  +p2[B12B2y+B3(1m)x](w1eδRt+w2eδLt)  +c2B2q2(w1eδRt+w2eδLt).
    $
    (7.22)

    Consider $ N_1 = (p_1q_1x-c_1)E_1 > 0 $ and $ N_2 = (p_2q_2y-c_2)E_2 > 0 $, we divide both sides of the first equation by that of the second equation of Eq (7.22)

    $ p1c1q1xp2c2q2y=B3(1m)y(p2c2q2y)+p1[A12A2KxA3(1m)y]+c1A2q1KA3(1m)x(p1c1q1x)+p2[B12B2y+B3(1m)x]+c2B2q2.
    $
    (7.23)

    Besides, the values of $ E_1 $ and $ E_2 $ at the interior equilibrium turn to

    $ E1=A1q1A2q1KxA3(1m)q1y,E2=B1q2B2q2y+B3(1m)q2x.
    $
    (7.24)

    Solving Eqs (7.18) (or $ (7.18)' $), (7.23) and (7.24) to get the optimal equilibrium solution $ x = x_{\tilde{\delta}}, y = y_{\tilde{\delta}} $ and the optimal harvesting efforts $ E_1 = E_{1\tilde{\delta}}, E_2 = E_{2\tilde{\delta}} $.

    In this section, some hypothetical data are taken to illustrate analytical results of the previous sections by mainly using the software MATLAB. Due to the features characterized by simulations should be analysed from qualitative perspective rather than quantitative one, then the parameters in system are not in view of real observed values. We consider three numerical examples to account for the stability of equilibria, the existence of bionomic equilibria as well as fuzzy optimal harvesting.

    Example 1 Consider the parameter values as follows: $ \tilde{r} = (2.00, 2.20, 2.40) $, $ \tilde{c} = (1.20, 1.35, 1.50) $, $ \tilde{d} = (0.30, 0.40, 0.50) $, $ \tilde{s} = (0.06, 0.07, 0.08) $, $ \tilde{e} = (0.60, 0.70, 0.80) $, $ K = 5 $ and $ m = 0.1 $.

    We discuss equilibria $ S_0 $, $ S_1 $ and $ S^* $ of the predator-prey model under different values of $ \alpha, w_1 $ and $ w_2 $, which are showed from Table 1 to Table 3.

    Table 1.  The trivial equilibrium $ S_0 $ for $ q_1 = 0.2, q_2 = 0.2, E_1 = 15 $ and $ E_2 = 10 $.
    $ w_1 $ $ w_2 $ $ S_0\ \mbox{at}\ \alpha=0 $ $ S_0\ \mbox{at}\ \alpha=0.3 $ $ S_0\ \mbox{at}\ \alpha=0.6 $ $ S_0\ \mbox{at}\ \alpha=0.9 $
    0 1 (0,0) (0,0) (0,0) (0,0)
    0.2 0.8 (0,0) (0,0) (0,0) (0,0)
    0.4 0.6 (0,0) (0,0) (0,0) (0,0)
    0.6 0.4 (0,0) (0,0) (0,0) (0,0)
    0.8 0.2 (0,0) (0,0) (0,0) (0,0)
    1 0 (0,0) (0,0) (0,0) (0,0)

     | Show Table
    DownLoad: CSV
    Table 2.  The axial equilibrium $ S_1 $ for $ q_1 = 0.2, q_2 = 0.2, E_1 = 5 $ and $ E_2 = 15 $.
    $ w_1 $ $ w_2 $ $ S_1\ \mbox{at}\ \alpha=0 $ $ S_1\ \mbox{at}\ \alpha=0.3 $ $ S_1\ \mbox{at}\ \alpha=0.6 $ $ S_1\ \mbox{at}\ \alpha=0.9 $
    0 1 (3.5000,0) (3.2524,0) (3.0189,0) (2.7982,0)
    0.2 0.8 (3.1731,0) (3.0340,0) (2.8996,0) (2.7697,0)
    0.4 0.6 (2.8704,0) (2.8269,0) (2.7839,0) (2.7413,0)
    0.6 0.4 (2.5893,0) (2.6302,0) (2.6715,0) (2.7132,0)
    0.8 0.2 (2.3276,0) (2.4431,0) (2.5623,0) (2.6854,0)
    1 0 (2.0833,0) (2.2650,0) (2.4561,0) (2.6577,0)

     | Show Table
    DownLoad: CSV
    Table 3.  The interior equilibrium $ S^* $ for $ q_1 = 0.2, q_2 = 0.2, E_1 = 5 $ and $ E_2 = 2 $.
    $ w_1 $ $ w_2 $ $ S^*\ \mbox{at}\ \alpha=0 $ $ S^*\ \mbox{at}\ \alpha=0.3 $ $ S^*\ \mbox{at}\ \alpha=0.6 $ $ S^*\ \mbox{at}\ \alpha=0.9 $
    0 1 (1.0479,0.9082) (1.1245,0.7824) (1.2067,0.6618) (1.2952,0.5454)
    0.2 0.8 (1.1513,0.7417) (1.2011,0.6697) (1.2531,0.5992) (1.3075,0.5301)
    0.4 0.6 (1.2650,0.5838) (1.2830,0.5607) (1.3013,0.5378) (1.3199,0.5150)
    0.6 0.4 (1.3907,0.4323) (1.3710,0.4547) (1.3516,0.4772) (1.3325,0.4998)
    0.8 0.2 (1.5307,0.2853) (1.4658,0.3510) (1.4041,0.4175) (1.3452,0.4847)
    1 0 (1.6875,0.1407) (1.5682,0.2490) (1.4588,0.3584) (1.3580,0.4697)

     | Show Table
    DownLoad: CSV

    Tables 13 demonstrate the trivial equilibrium is always fixed at $ (0, 0) $ with different combinations of $ w_1 $, $ w_2 $ and $ \alpha $, the values of prey $ x $ and predator $ y $ stay at zero; for the axial equilibrium, the value of prey $ x $ decreases under the same $ \alpha $ with increasing $ w_1 $, predator $ y $ is invariant in zero; for the interior equilibrium, the value of prey $ x $ increases under the same $ \alpha $ with increasing $ w_1 $, and predator $ y $ decreases under the same $ \alpha $ with increasing $ w_1 $.

    Figures 13 are time series of prey $ x $ and predator $ y $ with initial values $ (2, 6) $ for different $ w_1 $, $ w_2 $ and $ \alpha $, which show the same results as what Tables 1–3 reflect. The initial fluctuate for both species gradually saturate to a steady state level over time.

    Figure 1.  (a)–(f) Time series diagram of prey $ x $ and predator $ y $ with initial values $ (2, 6) $ for $ q_1 = 0.2, q_2 = 0.2, E_1 = 15, E_2 = 10 $ and different $ w_1, w_2 $ and $ \alpha $, respectively, $ t\in[0,100] $.
    Figure 2.  (a)–(f) Time series diagram of prey $ x $ and predator $ y $ with initial values $ (2, 6) $ for $ q_1 = 0.2, q_2 = 0.2, E_1 = 5, E_2 = 15 $ and different $ w_1, w_2 $ and $ \alpha $, respectively, $ t\in[0,100] $.
    Figure 3.  (a)–(f) Time series diagram of prey $ x $ and predator $ y $ with initial values $ (2, 6) $ for $ q_1 = 0.2, q_2 = 0.2, E_1 = 5, E_2 = 2 $ and different $ w_1, w_2 $ and $ \alpha $, respectively, $ t\in[0,100] $.

    The phase portrait of prey $ x $ and predator $ y $ corresponding to stabilities of different equilibria for all kinds of combinations of $ w_1 $, $ w_2 $ and $ \alpha $ are shown in Figures 46, respectively. Meanwhile, the equilibria are stable under different initial conditions.

    Figure 4.  (a)–(f) Phase portrait of prey $ x $ and predator $ y $ with initial values $ (2.5, 6), (3, 3) $ and $ (3, 5) $ for $ q_1 = 0.2, q_2 = 0.2, E_1 = 15, E_2 = 10 $ and different $ w_1 $, $ w_2 $ and $ \alpha $, respectively, $ t\in[0,100] $.
    Figure 5.  (a)–(f) Phase portrait of prey $ x $ and predator $ y $ with initial values $ (2.5, 6), (3, 3) $ and $ (3, 5) $ for $ q_1 = 0.2, q_2 = 0.2, E_1 = 5, E_2 = 15 $ and different $ w_1 $, $ w_2 $ and $ \alpha $, respectively, $ t\in[0,100] $.
    Figure 6.  (a)–(f) Phase portrait of prey $ x $ and predator $ y $ with initial values $ (2.5, 6), (2.5, 3) $ and $ (2.5, 5) $ for $ q_1 = 0.2, q_2 = 0.2, E_1 = 5, E_2 = 2 $ and different $ w_1 $, $ w_2 $ and $ \alpha $, respectively, $ t\in[0,100] $.

    Dynamical behavior of prey $ x $ and predator $ y $ under interior steady state level with respect to $ \alpha $ for initial condition $ (x(0), y(0)) = (2.5, 6) $, and $ q_1 = 0.2, q_2 = 0.2, E_1 = 5, E_2 = 2 $ are presented in Figure 7.

    Figure 7.  (a)–(f) Dynamical behavior of prey $ x $ and predator $ y $ with respect to $ \alpha $ for initial condition $ (x(0), y(0)) = (2.5, 6) $, and $ q_1 = 0.2, q_2 = 0.2, E_1 = 5, E_2 = 2 $.

    The figures in Figure 7 reflect the interior equilibrium changes for different $ \alpha $ which is in accordance with the results in Table 3. As $ \alpha $ increases, prey $ x $ increases and predator $ y $ decreases for $ w_1 = 0, w_2 = 1 $ and $ w_1 = 0.2, w_2 = 0.8 $ and $ w_1 = 0.4, w_2 = 0.6 $. However, in the cases $ w_1 = 0.6, w_2 = 0.4 $ and $ w_1 = 0.8, w_2 = 0.2 $ and $ w_1 = 1, w_2 = 0 $, as $ \alpha $ increases, prey $ x $ decreases and predator $ y $ increases. Therefore, the parameter $ \alpha $ plays an major part in the stability of interior equilibrium.

    In the mean time, the dynamical behavior of prey $ x $ and predator $ y $ about $ w_1 $ and $ w_2 $ under initial condition $ (x(0), y(0)) = (2.5, 6) $ for $ \alpha = 0, 0.3, 0.6 $ and $ \alpha = 0.9 $, and $ q_1 = 0.2, q_2 = 0.2, E_1 = 5, E_2 = 2 $ are presented in Figure 8. These figures depict the same results as that in Table 3, i.e., for the interior equilibrium, the value of prey $ x $ increases under fixed $ \alpha $ with increasing $ w_1 $, and predator $ y $ decreases under fixed $ \alpha $ with increasing $ w_1 $.

    Figure 8.  (a)–(d) Dynamical behavior of prey $ x $ and predator $ y $ with respect to $ w_1 $ and $ w_2 $ under initial condition $ (x(0), y(0)) = (2.5, 6) $ for $ \alpha = 0, 0.3, 0.6 $ and $ \alpha = 0.9 $ respectively, and $ q_1 = 0.2 $, $ q_2 = 0.2 $, $ E_1 = 5 $, $ E_2 = 2 $.

    Example 2 Consider the parameter values as follows: $ \tilde{r} = (1.50, 1.55, 1.60) $, $ \tilde{c} = (0.250, 0.275, 0.300) $, $ \tilde{d} = (0.450, 0.475, 0.500) $, $ \tilde{s} = (0.150, 0.175, 0.200) $, $ \tilde{e} = (1.300, 1.325, 1.350) $, $ K = 10, m = 0.1, q_1 = 0.92, q_2 = 0.95, p_1 = 20, p_2 = 25, c_1 = 30 $ and $ c_2 = 15 $.

    The bionomic equilibria for different combinations of $ w_1 $, $ w_2 $ and $ \alpha $ are emerged in Table 4. It is observed that $ x_\infty $ and $ y_\infty $ are invariable with respect to $ w_1 $, $ w_2 $ and $ \alpha $, and $ E_{1\infty} $ and $ E_{2\infty} $ are decreasing with fixed $ \alpha $ and increasing $ w_1 $. As well as, together with the value of $ \alpha $ increasing, $ E_{1\infty} $ and $ E_{2\infty} $ decreases for $ w_1 = 0, w_2 = 1 $ and $ w_1 = 0.2, w_2 = 0.8 $ and $ w_1 = 0.4, w_2 = 0.6 $. Nevertheless, in the cases $ w_1 = 0.6, w_2 = 0.4 $ and $ w_1 = 0.8, w_2 = 0.2 $ and $ w_1 = 1, w_2 = 0 $, as the value of $ \alpha $ increases, $ E_{1\infty} $ and $ E_{2\infty} $ are increasing.

    Table 4.  Positive bionomic equilibrium for different $ w_1, w_2 $ and $ \alpha $.
    $ w_1 $ $ w_2 $ $ \alpha=0 $ $ \alpha=0.3 $ $ \alpha=0.6 $ $ \alpha=0.9 $
    $ (x_\infty, y_\infty) $ $ (E_{1\infty}, E_{2\infty}) $ $ (x_\infty, y_\infty) $ $ (E_{1\infty}, E_{2\infty}) $ $ (x_\infty, y_\infty) $ $ (E_{1\infty}, E_{2\infty}) $ $ (x_\infty, y_\infty) $ $ (E_{1\infty}, E_{2\infty}) $
    0 1 (1.6304,0.6316) (1.3188,1.5118) (1.6304,0.6316) (1.2952,1.4874) (1.6304,0.6316) (1.2716,1.4629) (1.6304,0.6316) (1.2480,1.4384)
    0.2 0.8 (1.6304,0.6316) (1.2874,1.4792) (1.6304,0.6316) (1.2732,1.4645) (1.6304,0.6316) (1.2591,1.4499) (1.6304,0.6316) (1.2449,1.4352)
    0.4 0.6 (1.6304,0.6316) (1.2559,1.4466) (1.6304,0.6316) (1.2512,1.4417) (1.6304,0.6316) (1.2465,1.4368) (1.6304,0.6316) (1.2418,1.4319)
    0.6 0.4 (1.6304,0.6316) (1.2245,1.4140) (1.6304,0.6316) (1.2292,1.4189) (1.6304,0.6316) (1.2339,1.4238) (1.6304,0.6316) (1.2386,1.4287)
    0.8 0.2 (1.6304,0.6316) (1.1930,1.3814) (1.6304,0.6316) (1.2071,1.3960) (1.6304,0.6316) (1.2213,1.4107) (1.6304,0.6316) (1.2355,1.4254)
    1 0 (1.6304,0.6316) (1.1615,1.3487) (1.6304,0.6316) (1.1851,1.3732) (1.6304,0.6316) (1.2087,1.3977) (1.6304,0.6316) (1.2323,1.4221)

     | Show Table
    DownLoad: CSV

    Example 3 Consider the parameter values as follows: $ \tilde{r} = (2.80, 2.83, 2.85), \tilde{c} = (2.45, 2.55, 2.65), \tilde{d} = (0.012, 0.013, 0.015), \tilde{s} = (0.008, 0.009, 0.010), \tilde{e} = (0.20, 0.21, 0.22), K = 10, m = 0.1, q_1 = 0.45, q_2 = 0.55, p_1 = 30, p_2 = 25, c_1 = 25, c_2 = 12 $ and $ \tilde{\delta} = (0.8, 0.9, 1.0) $.

    The optimal equilibrium and harvesting effort for different combinations $ w_1, w_2, \alpha $ and $ w_1', w_2', \alpha' $ are appeared in Tables 59, respectively. For different combinations $ w_1, w_2, \alpha $ and $ w_1', w_2', \alpha' $, we see the existence of optimal equilibrium.

    Table 5.  Optimal equilibrium and optimal harvesting effort for $ w_1 = 0, w_2 = 1, \alpha = 0 $.
    $ w_1' $ $ w_2' $ $ \alpha'=0 $ $ \alpha'=0.3 $ $ \alpha'=0.6 $ $ \alpha'=0.9 $
    $ (x_{\tilde{\delta}}, y_{\tilde{\delta}}) $ $ (E_{1\tilde{\delta}}, E_{2\tilde{\delta}}) $ $ (x_{\tilde{\delta}}, y_{\tilde{\delta}}) $ $ (E_{1\tilde{\delta}}, E_{2\tilde{\delta}}) $ $ (x_{\tilde{\delta}}, y_{\tilde{\delta}}) $ $ (E_{1\tilde{\delta}}, E_{2\tilde{\delta}}) $ $ (x_{\tilde{\delta}}, y_{\tilde{\delta}}) $ $ (E_{1\tilde{\delta}}, E_{2\tilde{\delta}}) $
    0.2 0.8 (1.8964,1.0333) (0.0902,0.6461) (1.8966,1.0332) (0.0906,0.6460) (1.8966,1.0331) (0.0908,0.6460) (1.8967,1.0331) (0.0909,0.6460)
    0.4 0.6 (1.8964,1.0333) (0.0900,0.6463) (1.8965,1.0332) (0.0904,0.6461) (1.8966,1.0331) (0.0907,0.6460) (1.8967,1.0331) (0.0909,0.6460)
    0.5 0.5 (1.8964,1.0333) (0.0900,0.6463) (1.8965,1.0332) (0.0904,0.6461) (1.8967,1.0331) (0.0908,0.6460) (1.8967,1.0331) (0.0908,0.6460)
    0.6 0.4 (1.8964,1.0333) (0.0900,0.6463) (1.8965,1.0332) (0.0904,0.6461) (1.8966,1.0331) (0.0907,0.6460) (1.8967,1.0331) (0.0908,0.6460)
    0.8 0.2 (1.8964,1.0333) (0.0902,0.6461) (1.8966,1.0332) (0.0906,0.6460) (1.8966,1.0331) (0.0908,0.6460) (1.8967,1.0331) (0.0909,0.6460)

     | Show Table
    DownLoad: CSV
    Table 6.  Optimal equilibrium and optimal harvesting effort for $ w_1 = 0, w_2 = 1, \alpha = 0.5 $.
    $ w_1' $ $ w_2' $ $ \alpha'=0 $ $ \alpha'=0.3 $ $ \alpha'=0.6 $ $ \alpha'=0.9 $
    $ (x_{\tilde{\delta}}, y_{\tilde{\delta}}) $ $ (E_{1\tilde{\delta}}, E_{2\tilde{\delta}}) $ $ (x_{\tilde{\delta}}, y_{\tilde{\delta}}) $ $ (E_{1\tilde{\delta}}, E_{2\tilde{\delta}}) $ $ (x_{\tilde{\delta}}, y_{\tilde{\delta}}) $ $ (E_{1\tilde{\delta}}, E_{2\tilde{\delta}}) $ $ (x_{\tilde{\delta}}, y_{\tilde{\delta}}) $ $ (E_{1\tilde{\delta}}, E_{2\tilde{\delta}}) $
    0.2 0.8 (1.9230,0.9964) (0.1262,0.6384) (1.9232,0.9963) (0.1264,0.6385) (1.9234,0.9963) (0.1265,0.6386) (1.9235,0.9963) (0.1266,0.6386)
    0.4 0.6 (1.9230,0.9964) (0.1262,0.6384) (1.9232,0.9963) (0.1264,0.6385) (1.9234,0.9963) (0.1265,0.6386) (1.9235,0.9963) (0.1266,0.6386)
    0.5 0.5 (1.9230,0.9964) (0.1262,0.6384) (1.9232,0.9963) (0.1264,0.6385) (1.9234,0.9963) (0.1265,0.6386) (1.9235,0.9963) (0.1266,0.6386)
    0.6 0.4 (1.9230,0.9964) (0.1262,0.6384) (1.9232,0.9963) (0.1264,0.6385) (1.9234,0.9963) (0.1265,0.6386) (1.9235,0.9963) (0.1266,0.6386)
    0.8 0.2 (1.9230,0.9964) (0.1262,0.6384) (1.9232,0.9963) (0.1264,0.6385) (1.9234,0.9963) (0.1265,0.6386) (1.9235,0.9963) (0.1266,0.6386)

     | Show Table
    DownLoad: CSV
    Table 7.  Optimal equilibrium and optimal harvesting effort for $ w_1 = 0, w_2 = 1, \alpha = 0.9 $.
    $ w_1' $ $ w_2' $ $ \alpha'=0 $ $ \alpha'=0.3 $ $ \alpha'=0.6 $ $ \alpha'=0.9 $
    $ (x_{\tilde{\delta}}, y_{\tilde{\delta}}) $ $ (E_{1\tilde{\delta}}, E_{2\tilde{\delta}}) $ $ (x_{\tilde{\delta}}, y_{\tilde{\delta}}) $ $ (E_{1\tilde{\delta}}, E_{2\tilde{\delta}}) $ $ (x_{\tilde{\delta}}, y_{\tilde{\delta}}) $ $ (E_{1\tilde{\delta}}, E_{2\tilde{\delta}}) $ $ (x_{\tilde{\delta}}, y_{\tilde{\delta}}) $ $ (E_{1\tilde{\delta}}, E_{2\tilde{\delta}}) $
    0.2 0.8 (1.9361,0.9744) (0.1272,0.6293) (1.9366,0.9744) (0.1270,0.6294) (1.9369,0.9744) (0.1268,0.6295) (1.9370,0.9744) (0.1267,0.6296)
    0.4 0.6 (1.9361,0.9744) (0.1272,0.6293) (1.9366,0.9744) (0.1270,0.6294) (1.9369,0.9744) (0.1268,0.6295) (1.9370,0.9744) (0.1267,0.6296)
    0.5 0.5 (1.9361,0.9744) (0.1272,0.6293) (1.9366,0.9744) (0.1270,0.6294) (1.9369,0.9744) (0.1268,0.6295) (1.9370,0.9744) (0.1267,0.6296)
    0.6 0.4 (1.9361,0.9744) (0.1272,0.6293) (1.9366,0.9744) (0.1270,0.6294) (1.9369,0.9744) (0.1268,0.6295) (1.9370,0.9744) (0.1267,0.6296)
    0.8 0.2 (1.9361,0.9744) (0.1272,0.6293) (1.9366,0.9744) (0.1270,0.6294) (1.9369,0.9744) (0.1268,0.6295) (1.9370,0.9744) (0.1267,0.6296)

     | Show Table
    DownLoad: CSV
    Table 8.  Optimal equilibrium and optimal harvesting effort for $ w_1 = 0.5, w_2 = 0.5, \alpha = 0.9 $.
    $ w_1' $ $ w_2' $ $ \alpha'=0 $ $ \alpha'=0.3 $ $ \alpha'=0.6 $ $ \alpha'=0.9 $
    $ (x_{\tilde{\delta}}, y_{\tilde{\delta}}) $ $ (E_{1\tilde{\delta}}, E_{2\tilde{\delta}}) $ $ (x_{\tilde{\delta}}, y_{\tilde{\delta}}) $ $ (E_{1\tilde{\delta}}, E_{2\tilde{\delta}}) $ $ (x_{\tilde{\delta}}, y_{\tilde{\delta}}) $ $ (E_{1\tilde{\delta}}, E_{2\tilde{\delta}}) $ $ (x_{\tilde{\delta}}, y_{\tilde{\delta}}) $ $ (E_{1\tilde{\delta}}, E_{2\tilde{\delta}}) $
    0.2 0.8 (1.9406,0.9694) (0.1236,0.6273) (1.9417,0.9695) (0.1224,0.6277) (1.9429,0.9696) (0.1212,0.6280) (1.9440,0.9697) (0.1198,0.6284)
    0.4 0.6 (1.9406,0.9694) (0.1236,0.6273) (1.9417,0.9695) (0.1224,0.6277) (1.9429,0.9696) (0.1212,0.6280) (1.9440,0.9697) (0.1198,0.6284)
    0.5 0.5 (1.9406,0.9694) (0.1236,0.6273) (1.9417,0.9695) (0.1224,0.6277) (1.9429,0.9696) (0.1212,0.6280) (1.9440,0.9697) (0.1198,0.6284)
    0.6 0.4 (1.9406,0.9694) (0.1236,0.6273) (1.9417,0.9695) (0.1224,0.6277) (1.9429,0.9696) (0.1212,0.6280) (1.9440,0.9697) (0.1198,0.6284)
    0.8 0.2 (1.9406,0.9694) (0.1236,0.6273) (1.9417,0.9695) (0.1224,0.6277) (1.9429,0.9696) (0.1212,0.6280) (1.9440,0.9697) (0.1198,0.6284)

     | Show Table
    DownLoad: CSV
    Table 9.  Optimal equilibrium and optimal harvesting effort for $ w_1 = 1, w_2 = 0, \alpha = 0.9 $.
    $ w_1' $ $ w_2' $ $ \alpha'=0 $ $ \alpha'=0.3 $ $ \alpha'=0.6 $ $ \alpha'=0.9 $
    $ (x_{\tilde{\delta}}, y_{\tilde{\delta}}) $ $ (E_{1\tilde{\delta}}, E_{2\tilde{\delta}}) $ $ (x_{\tilde{\delta}}, y_{\tilde{\delta}}) $ $ (E_{1\tilde{\delta}}, E_{2\tilde{\delta}}) $ $ (x_{\tilde{\delta}}, y_{\tilde{\delta}}) $ $ (E_{1\tilde{\delta}}, E_{2\tilde{\delta}}) $ $ (x_{\tilde{\delta}}, y_{\tilde{\delta}}) $ $ (E_{1\tilde{\delta}}, E_{2\tilde{\delta}}) $
    0.2 0.8 (1.9435,0.9648) (0.1196,0.6246) (1.9434,0.9648) (0.1197,0.6246) (1.9433,0.9648) (0.1197,0.6246) (1.9433,0.9648) (0.1197,0.6247)
    0.4 0.6 (1.9436,0.9648) (0.1197,0.6245) (1.9434,0.9648) (0.1197,0.6246) (1.9434,0.9648) (0.1197,0.6246) (1.9433,0.9648) (0.1197,0.6246)
    0.5 0.5 (1.9436,0.9648) (0.1197,0.6245) (1.9435,0.9648) (0.1197,0.6246) (1.9434,0.9648) (0.1197,0.6246) (1.9433,0.9648) (0.1197,0.6246)
    0.6 0.4 (1.9436,0.9648) (0.1197,0.6245) (1.9434,0.9648) (0.1197,0.6246) (1.9434,0.9648) (0.1197,0.6246) (1.9433,0.9648) (0.1197,0.6247)
    0.8 0.2 (1.9435,0.9648) (0.1196,0.6246) (1.9434,0.9648) (0.1197,0.6246) (1.9433,0.9648) (0.1197,0.6246) (1.9433,0.9648) (0.1197,0.6247)

     | Show Table
    DownLoad: CSV

    In Table 5, when $ \alpha' = 0 $ it is observed that prey $ x_{\tilde{\delta}} $ and predator $ y_{\tilde{\delta}} $ are invariant, optimal harvesting effort of the prey $ E_{1\tilde{\delta}} $ decreases and then increases, but the optimal harvesting effort of the predator $ E_{2\tilde{\delta}} $ increases and then decreases with increasing $ w_1' $; when $ \alpha' = 0.3 $, $ x_{\tilde{\delta}} $ decreases and then increases, $ y_{\tilde{\delta}} $ is invariant, $ E_{1\tilde{\delta}} $ decreases and then increases, and $ E_{2\tilde{\delta}} $ increases and then decreases with increasing $ w_1' $; when $ \alpha' = 0.6 $, $ x_{\tilde{\delta}} $ increases and then decreases, $ y_{\tilde{\delta}} $ is invariant, $ E_{1\tilde{\delta}} $ fluctuates between 0.0907 and 0.0908, and $ E_{2\tilde{\delta}} $ is invariant with increasing $ w_1' $; when $ \alpha' = 0.9 $, $ x_{\tilde{\delta}} $ and $ y_{\tilde{\delta}} $ are invariant, $ E_{1\tilde{\delta}} $ decreases and then increases, and $ E_{2\tilde{\delta}} $ is invariant with increasing $ w_1' $. Besides, when fixed $ w_1' $ and $ w_2' $, $ x_{\tilde{\delta}} $ increases, $ y_{\tilde{\delta}} $ decreases, $ E_{1\tilde{\delta}} $ increases, and $ E_{2\tilde{\delta}} $ decreases with the increase of $ \alpha' $.

    In Table 6, when fixed $ \alpha' $, $ x_{\tilde{\delta}} $, $ y_{\tilde{\delta}} $, $ E_{1\tilde{\delta}} $ and $ E_{2\tilde{\delta}} $ are invariant with the increase of $ w_1' $; when fixed $ w_1' $ and $ w_2' $, $ x_{\tilde{\delta}} $ increases, $ y_{\tilde{\delta}} $ decreases, $ E_{1\tilde{\delta}} $ increases, and $ E_{2\tilde{\delta}} $ increases with the increase of $ \alpha' $.

    In Table 7, when fixed $ \alpha' $, $ x_{\tilde{\delta}} $, $ y_{\tilde{\delta}} $, $ E_{1\tilde{\delta}} $ and $ E_{2\tilde{\delta}} $ are invariant with the increase of $ w_1' $; when fixed $ w_1' $ and $ w_2' $, $ x_{\tilde{\delta}} $ increases, $ y_{\tilde{\delta}} $ is invariant, $ E_{1\tilde{\delta}} $ decreases, and $ E_{2\tilde{\delta}} $ increases with the increase of $ \alpha' $.

    In Table 8, when fixed $ \alpha' $, $ x_{\tilde{\delta}} $, $ y_{\tilde{\delta}} $, $ E_{1\tilde{\delta}} $ and $ E_{2\tilde{\delta}} $ are invariant with the increase of $ w_1' $; when fixed $ w_1' $ and $ w_2' $, $ x_{\tilde{\delta}} $ increases, $ y_{\tilde{\delta}} $ increases, $ E_{1\tilde{\delta}} $ decreases, and $ E_{2\tilde{\delta}} $ increases with the increase of $ \alpha' $.

    In Table 9, when $ \alpha' = 0 $, $ x_{\tilde{\delta}} $ increases and then decreases and $ y_{\tilde{\delta}} $ is invariant, $ E_{1\tilde{\delta}} $ increases and then decreases, but $ E_{2\tilde{\delta}} $ decreases and then increases with increasing $ w_1' $; when $ \alpha' = 0.3 $ or $ \alpha' = 0.6 $, we can see that $ x_{\tilde{\delta}} $ increases and then decreases, $ y_{\tilde{\delta}} $, $ E_{1\tilde{\delta}} $ and $ E_{2\tilde{\delta}} $ are invariant with increasing $ w_1' $; when $ \alpha' = 0.9 $, $ x_{\tilde{\delta}} $, $ y_{\tilde{\delta}} $ and $ E_{1\tilde{\delta}} $ are invariant, $ E_{2\tilde{\delta}} $ decreases and then increases with increasing $ w_1' $. When fixed $ w_1' $ and $ w_2' $, $ x_{\tilde{\delta}} $ decreases, $ y_{\tilde{\delta}} $ is invariant, $ E_{1\tilde{\delta}} $ is almost invariant, and $ E_{2\tilde{\delta}} $ increases with the increase of $ \alpha' $.

    From Tables 57, it is easy to see that when fixed $ w_1, w_2 $ and $ w_1', w_2', \alpha' $, the prey increases, the predator decreases, the optimal harvesting effort of the prey increases and the optimal harvesting effort of the predator decreases with increasing $ \alpha $.

    Based on Tables 79, we observed that when fixed $ \alpha $ and $ w_1', w_2', \alpha' $, the prey increases, the predator decreases, the optimal harvesting effort of the prey decreases and the optimal harvesting effort of the predator decreases with increasing $ w_1 $.

    In the field of biomathematics most harvesting models are on account of the assumption with precise parameters, but any species in nature will inevitably be affected by the complexity of the ecosystem itself and the limitations of people's cognition leading to the uncertainty of the biological parameters. These uncertainties can be roughly classified into three categories: fuzziness, randomness, and unascertainty. In view of this point, we have modified a fuzzy parameters predator-prey model with refuge. Here, the fuzzy parameters are introduced into biotic potential of the prey ($ r $), predator mortality ($ d $), intraspecific competition rate of the predator ($ s $), predation rate ($ c $) and the product of predation rate and conversion rate ($ e $).

    We have discussed the positivity and boundedness of system (Theorem 4.1) and investigated the existence and stability of equilibria (Theorems 5.1 and 5.2). The tables and graphs for different equilibria are showed from Table 1 to Table 3 and from Figure 1 to Figure 3, respectively. Figures 46 reflect phase portrait of prey and predator about different $ w_1, w_2 $ and $ \alpha $. Dynamical behaviors of prey and predator about $ \alpha $ or $ w_1 $ and $ w_2 $, separately, are presented in Figures 7 and 8. We have also analysed the existence of bionomic equilibria (Theorem 6.1) and the mathematical results for positive bionomic equilibrium are displayed in Table 4.

    The highlight of this article is to solve the optimal control problem with fuzzy inflation and discount, that is maximize the fuzzy objective functional $ \tilde{J} $. The optimal steady state solution is displayed in Tables 59. We obtain optimal equilibrium and optimal harvesting effort as respect to $ w_1', w_2' $ and $ \alpha' $ under five cases including $ w_1 = 0 $, $ w_2 = 1 $, $ \alpha = 0 $ and $ w_1 = 0 $, $ w_2 = 1 $, $ \alpha = 0.5 $ and $ w_1 = 0 $, $ w_2 = 1 $, $ \alpha = 0.9 $ and $ w_1 = 0.5 $, $ w_2 = 0.5 $, $ \alpha = 0.9 $ and $ w_1 = 1 $, $ w_2 = 0 $, $ \alpha = 0.9 $.

    In this paper, we consider the fuzzification of biological parameters in model. However, the fuzzification of unit capture cost or unit biomass price for both species seems to be also feasible. Future work will focus on the imprecision of economic parameters to enhance our model. Finally, this paper uses $ \alpha $-cut of triangular fuzzy number to describe imprecise parameters in model, however other common forms of fuzzy number such as normal fuzzy number and trapezoidal fuzzy number are important in fuzzy set theory, and then we are going to apply other forms of fuzzy number depicting imprecise parameters in model.

    The authors thank the editor and referees for their careful reading and valuable comments. The work is supported by the National Natural Science Foundation of China (No.12101211) and the Program for Innovative Research Team of the Higher Education Institution of Hubei Province (No.T201812) and the Scientific Research Project of Education Department of Hubei Province (No.B2020090).

    The authors have no conflict of interest.

    [1] Geiman TM, Robertson KD (2002) Chromatin remodeling, histone modifications, and DNA methylation-how does it all fit together? J Cell Biochem 87: 117–125. doi: 10.1002/jcb.10286
    [2] Bonev B, Cavalli G (2016) Organization and function of the 3D genome. Nat Rev Genet 17: 661–678. doi: 10.1038/nrg.2016.112
    [3] Bannister AJ, Kouzarides T (2011) Regulation of chromatin by histone modifications. Cell Res 21: 381–395. doi: 10.1038/cr.2011.22
    [4] Tessarz P, Kouzarides T (2014) Histone core modifications regulating nucleosome structure and dynamics. Nat Rev Mol Cell Bio 15: 703–708.
    [5] Koch A, Joosten SC, Feng Z, et al. (2018) Analysis of DNA methylation in cancer: Location revisited. Nat Rev Clin Oncol 15: 459–466. doi: 10.1038/s41571-018-0004-4
    [6] Watt F, Molloy PL (1988) Cytosine methylation prevents binding to DNA of a HeLa cell transcription factor required for optimal expression of the adenovirus major late promoter. Gene Dev 2: 1136–1143. doi: 10.1101/gad.2.9.1136
    [7] Csankovszki G, Nagy A, Jaenisch R (2001) Synergism of Xist RNA, DNA methylation, and histone hypoacetylation in maintaining X chromosome inactivation. J Cell Biol 153: 773–784. doi: 10.1083/jcb.153.4.773
    [8] Li E, Beard C, Jaenisch R (1993) Role for DNA methylation in genomic imprinting. Nature 366: 362–365. doi: 10.1038/366362a0
    [9] Okano M, Bell DW, Haber DA, et al. (1999) DNA methyltransferases Dnmt3a and Dnmt3b are essential for de novo methylation and mammalian development. Cell 99: 247–257. doi: 10.1016/S0092-8674(00)81656-6
    [10] Hashimoto H, Liu Y, Upadhyay AK, et al. (2012) Recognition and potential mechanisms for replication and erasure of cytosine hydroxymethylation. Nucleic Acids Res 40: 4841–4849. doi: 10.1093/nar/gks155
    [11] Nabel CS, Jia H, Ye Y, et al. (2012) AID/APOBEC deaminases disfavor modified cytosines implicated in DNA demethylation. Nat Chem Biol 8: 751–758. doi: 10.1038/nchembio.1042
    [12] Bochtler M, Kolano A, Xu GL (2017) DNA demethylation pathways: Additional players and regulators. Bioessays 39: 1–13. doi: 10.1002/bies.201670013
    [13] Tahiliani M, Koh KP, Shen Y, et al. (2009) Conversion of 5-methylcytosine to 5-hydroxymethylcytosine in mammalian DNA by MLL partner TET1. Science 324: 930–935. doi: 10.1126/science.1170116
    [14] Tamanaha E, Guan S, Marks K, et al. (2016) Distributive processing by the iron(II)/α-ketoglutarate-dependent catalytic domains of the TET enzymes is consistent with epigenetic roles for oxidized 5-methylcytosine bases. J Am Chem Soc 138: 9345–9348. doi: 10.1021/jacs.6b03243
    [15] Müller U, Bauer C, Siegl M, et al. (2014) TET-mediated oxidation of methylcytosine causes TDG or NEIL glycosylase dependent gene reactivation. Nucleic Acids Res 42: 8592–8604. doi: 10.1093/nar/gku552
    [16] Weber AR, Krawczyk C, Robertson AB, et al. (2016) Biochemical reconstitution of TET1-TDG-BER-dependent active DNA demethylation reveals a highly coordinated mechanism. Nat Commun 7: 10806. doi: 10.1038/ncomms10806
    [17] Wang J, Hevi S, Kurash JK, et al. (2009) The lysine demethylase LSD1 (KDM1) is required for maintenance of global DNA methylation. Nat Genet 41: 125–129. doi: 10.1038/ng.268
    [18] Szwagierczak A, Bultmann S, Schmidt CS, et al. (2010) Sensitive enzymatic quantification of 5-hydroxymethylcytosine in genomic DNA. Nucleic Acids Res 38: e181. doi: 10.1093/nar/gkq684
    [19] Spruijt CG, Gnerlich F, Smits AH, et al. (2013) Dynamic readers for 5-(hydroxy)methylcytosine and its oxidized derivatives. Cell 152: 1146–1159. doi: 10.1016/j.cell.2013.02.004
    [20] Matarese F, Carrillode SPE, Stunnenberg HG (2011) 5-Hydroxymethylcytosine: A new kid on the epigenetic block? Mol Syst Biol 7: 562.
    [21] Koh KP, Yabuuchi A, Rao S, et al. (2011) Tet1 and Tet2 regulate 5-hydroxymethylcytosine production and cell lineage specification in mouse embryonic stem cells. Cell Stem Cell 8: 200–213. doi: 10.1016/j.stem.2011.01.008
    [22] Ito S, D'Alessio AC, Taranova OV, et al. (2010) Role of Tet proteins in 5mC to 5hmC conversion, ES-cell self-renewal and inner cell mass specification. Nature 466: 1129–1133. doi: 10.1038/nature09303
    [23] Caron G, Hussein M, Kulis M, et al. (2015) Cell-cycle-dependent reconfiguration of the DNA methylome during terminal differentiation of human B cells into plasma cells. Cell Rep 13: 1059–1071. doi: 10.1016/j.celrep.2015.09.051
    [24] Costa Y, Ding J, Theunissen TW, et al. (2013) NANOG-dependent function of TET1 and TET2 in establishment of pluripotency. Nature 495: 370–374. doi: 10.1038/nature11925
    [25] Kriaucionis S, Heintz N (2009) The nuclear DNA base 5-hydroxymethylcytosine is present in Purkinje neurons and the brain. Science 324: 929–930. doi: 10.1126/science.1169786
    [26] Jin SG, Wu X, Li AX, et al. (2011) Genomic mapping of 5-hydroxymethylcytosine in the human brain. Nucleic Acids Res 39: 5015–5024. doi: 10.1093/nar/gkr120
    [27] Nestor CE, Ottaviano R, Reddington J, et al. (2012) Tissue type is a major modifier of the 5-hydroxymethylcytosine content of human genes. Genome Res 22: 467–477. doi: 10.1101/gr.126417.111
    [28] Szulwach KE, Li X, Li Y, et al. (2011) 5-hmC-mediated epigenetic dynamics during postnatal neurodevelopment and aging. Nat Neurosci 14: 1607–1616. doi: 10.1038/nn.2959
    [29] Sherwani SI, Khan HA (2015) Role of 5-hydroxymethylcytosine in neurodegeneration. Gene 570: 17–24. doi: 10.1016/j.gene.2015.06.052
    [30] Jeschke J, Collignon E, Fuks F (2016) Portraits of TET-mediated DNA hydroxymethylation in cancer. Curr Opin Genet Dev 36: 16–26. doi: 10.1016/j.gde.2016.01.004
    [31] Haffner MC, Chaux A, Meeker AK, et al. (2011) Global 5-hydroxymethylcytosine content is significantly reduced in tissue stem/progenitor cell compartments and in human cancers. Oncotarget 2: 627–637.
    [32] Mahé EA, Madigou T, Sérandour AA, et al. (2017) Cytosine modifications modulate the chromatin architecture of transcriptional enhancers. Genome Res 27: 947–958. doi: 10.1101/gr.211466.116
    [33] Nan X, Ng HH, Johnson CA, et al. (1998) Transcriptional repression by the methyl-CpG-binding protein MeCP2 involves a histone deacetylase complex. Nature 393: 386–389. doi: 10.1038/30764
    [34] Ng HH, Bird A (1999) DNA methylation and chromatin modification. Curr Opin Genet Dev 9: 158–163. doi: 10.1016/S0959-437X(99)80024-0
    [35] Schübeler D (2015) Function and information content of DNA methylation. Nature 517: 321–326. doi: 10.1038/nature14192
    [36] Williams K, Christensen J, Pedersen MT, et al. (2011) TET1 and hydroxymethylcytosine in transcription and DNA methylation fidelity. Nature 473: 343–348. doi: 10.1038/nature10066
    [37] Sérandour AA, Avner S, Oger F, et al. (2012) Dynamic hydroxymethylation of deoxyribonucleic acid marks differentiation-associated enhancers. Nucleic Acids Res 40: 8255–8265. doi: 10.1093/nar/gks595
    [38] Song CX, He C (2013) Potential functional roles of DNA demethylation intermediates. Trends Biochem Sci 38: 480–484. doi: 10.1016/j.tibs.2013.07.003
    [39] Sepulveda H, Villagra A, Montecino M (2017) Tet-mediated DNA demethylation is required for SWI/SNF-dependent chromatin remodeling and histone-modifying activities that trigger expression of the Sp7 osteoblast master gene during mesenchymal lineage commitment. Mol Cell Biol 37.
    [40] Yildirim O, Li R, Hung JH, et al. (2011) Mbd3/NURD complex regulates expression of 5-hydroxymethylcytosine marked genes in embryonic stem cells. Cell 147: 1498–1510. doi: 10.1016/j.cell.2011.11.054
    [41] Neri F, Incarnato D, Krepelova A, et al. (2013) Genome-wide analysis identifies a functional association of Tet1 and Polycomb repressive complex 2 in mouse embryonic stem cells. Genome Biol 14: R91. doi: 10.1186/gb-2013-14-8-r91
    [42] Deplus R, Delatte B, Schwinn MK, et al. (2013) TET2 and TET3 regulate GlcNAcylation and H3K4 methylation through OGT and SET1/COMPASS. EMBO J 32: 645–655. doi: 10.1038/emboj.2012.357
    [43] Kong L, Tan L, Lv R, et al. (2016) A primary role of TET proteins in establishment and maintenance of De Novo bivalency at CpG islands. Nucleic Acids Res 44: 8682–8692. doi: 10.1093/nar/gkw529
    [44] Mendonca A, Chang EH, Liu W, et al. (2014) Hydroxymethylation of DNA influences nucleosomal conformation and stability in vitro. BBA-Gene Regul Mech 1839: 1323–1329.
    [45] Deplus R, Delatte B, Schwinn MK, et al. (2013) TET2 and TET3 regulate GlcNAcylation and H3K4 methylation through OGT and SET1/COMPASS. EMBO J 32: 645–655. doi: 10.1038/emboj.2012.357
    [46] Guan W, Guyot R, Samarut J, et al. (2017) Methylcytosine dioxygenase TET3 interacts with thyroid hormone nuclear receptors and stabilizes their association to chromatin. P Natl Acad Sci USA 114: 8229–8234. doi: 10.1073/pnas.1702192114
    [47] Zhang YW, Wang Z, Xie W, et al. (2017) Acetylation enhances TET2 function in protecting against abnormal DNA methylation during oxidative stress. Mol Cell 65: 323–335. doi: 10.1016/j.molcel.2016.12.013
    [48] Cimmino L, Dawlaty MM, Ndiaye-Lobry D, et al. (2015) TET1 is a tumor suppressor of hematopoietic malignancy. Nat Immunol 16: 653–662. doi: 10.1038/ni.3148
    [49] An J, González-Avalos E, Chawla A, et al. (2015) Acute loss of TET function results in aggressive myeloid cancer in mice. Nat Commun 6: 10071. doi: 10.1038/ncomms10071
    [50] Kafer GR, Li X, Horii T, et al. (2016) 5-Hydroxymethylcytosine marks sites of DNA damage and promotes genome stability. Cell Rep 14: 1283–1292. doi: 10.1016/j.celrep.2016.01.035
    [51] Mahfoudhi E, Talhaoui I, Cabagnols X, et al. (2016) TET2-mediated 5-hydroxymethylcytosine induces genetic instability and mutagenesis. DNA Rep 43: 78–88. doi: 10.1016/j.dnarep.2016.05.031
    [52] Jiang D, Zhang Y, Hart RP, et al. (2015) Alteration in 5-hydroxymethylcytosine-mediated epigenetic regulation leads to Purkinje cell vulnerability in ATM deficiency. Brain 138: 3520–3536. doi: 10.1093/brain/awv284
    [53] Jiang D, Wei S, Chen F, et al. (2017) TET3-mediated DNA oxidation promotes ATR-dependent DNA damage response. EMBO Rep 18: 781–796. doi: 10.15252/embr.201643179
    [54] Blackford AN, Jackson SP (2017) ATM, ATR, and DNA-PK: the trinity at the heart of the DNA damage response. Mol Cell 66: 801–817. doi: 10.1016/j.molcel.2017.05.015
    [55] Sellou H, Lebeaupin T, Chapuis C, et al. (2016) The poly(ADP-ribose)-dependent chromatin remodeler Alc1 induces local chromatin relaxation upon DNA damage. Mol Biol Cell 27: 3791–3799. doi: 10.1091/mbc.e16-05-0269
    [56] Smith R, Sellou H, Chapuis C, et al. (2018) CHD3 and CHD4 recruitment and chromatin remodeling activity at DNA breaks is promoted by early poly(ADP-ribose)-dependent chromatin relaxation. Nucleic Acids Res 46: 6087.
    [57] Ciccarone F, Valentini E, Zampieri M, et al. (2015) 5mC-hydroxylase activity is influenced by the PARylation of TET1 enzyme. Oncotarget 6: 24333–24347.
    [58] Ciccarone F, Valentini E, Bacalini MG, et al. (2014) Poly(ADP-ribosyl)ation is involved in the epigenetic control of TET1 gene transcription. Oncotarget 5: 10356–10367.
    [59] Chou DM, Adamson B, Dephoure NE, et al. (2010) A chromatin localization screen reveals poly (ADP ribose)-regulated recruitment of the repressive polycomb and NuRD complexes to sites of DNA damage. P Natl Acad Sci USA 107: 18475–18480. doi: 10.1073/pnas.1012946107
    [60] Luijsterburg MS, Dinant C, Lans H, et al. (2009) Heterochromatin protein 1 is recruited to various types of DNA damage. J Cell Biol 185: 577–586. doi: 10.1083/jcb.200810035
    [61] Abu-Zhayia ER, Awwad SW, Ben-Oz B, et al. (2017) CDYL1 fosters double-strand break-induced transcription silencing and promotes homology-directed repair. J Mol Cell Biol 1: 1.
    [62] D'Alessandro G, Fagagna FDD (2017) Transcription and DNA damage: holding hands or crossing swords? J Mol Biol 429: 3215–3229. doi: 10.1016/j.jmb.2016.11.002
    [63] Puc J, Aggarwal AK, Rosenfeld MG (2017) Physiological functions of programmed DNA breaks in signal-induced transcription. Nat Rev Mol Cell Bio 18: 471–476.
    [64] Ju BG, Lunyak VV, Perissi V, et al. (2006) A topoisomerase IIbeta-mediated dsDNA break required for regulated transcription. Science 312: 1798–1802. doi: 10.1126/science.1127196
    [65] Madabhushi R, Gao F, Pfenning AR, et al. (2015) Activity-induced DNA breaks govern the expression of neuronal early-response genes. Cell 161: 1592–1605. doi: 10.1016/j.cell.2015.05.032
    [66] Perillo B, Ombra MN, Bertoni A, et al. (2008) DNA oxidation as triggered by H3K9me2 demethylation drives estrogen-induced gene expression. Science 319: 202–206. doi: 10.1126/science.1147674
    [67] Puc J, Kozbial P, Li W, et al. (2015) Ligand-dependent enhancer activation regulated by topoisomerase-I activity. Cell 160: 367–380. doi: 10.1016/j.cell.2014.12.023
    [68] Baranello L, Wojtowicz D, Cui K, et al. (2016) RNA polymerase II regulates topoisomerase 1 activity to favor efficient transcription. Cell 165: 357–371. doi: 10.1016/j.cell.2016.02.036
    [69] Bunch H, Lawney BP, Lin YF, et al. (2015) Transcriptional elongation requires DNA break-induced signalling. Nat Commun 6: 10191. doi: 10.1038/ncomms10191
    [70] Marnef A, Cohen S, Legube G (2017) Transcription-coupled DNA double-strand break repair: active genes need special care. J Mol Biol 429: 1277–1288. doi: 10.1016/j.jmb.2017.03.024
    [71] Huertas P, Aguilera A (2003) Cotranscriptionally formed DNA:RNA hybrids mediate transcription elongation impairment and transcription-associated recombination. Mol Cell 12: 711–721. doi: 10.1016/j.molcel.2003.08.010
    [72] Sollier J, Stork CT, García-Rubio ML, et al. (2014) Transcription-coupled nucleotide excision repair factors promote R-loop-induced genome instability. Mol Cell 56: 777–785. doi: 10.1016/j.molcel.2014.10.020
    [73] Métivier R, Gallais R, Tiffoche C, et al. (2008) Cyclical DNA methylation of a transcriptionally active promoter. Nature 452: 45–50. doi: 10.1038/nature06544
    [74] Li J, Wu X, Zhou Y, et al. (2018) Decoding the dynamic DNA methylation and hydroxymethylation landscapes in endodermal lineage intermediates during pancreatic differentiation of hESC. Nucleic Acids Res 46: 2883–2900. doi: 10.1093/nar/gky063
    [75] Zhang Y, Zhang D, Li Q, et al. (2016) Nucleation of DNA repair factors by FOXA1 links DNA demethylation to transcriptional pioneering. Nat Genet 48: 1003–1013. doi: 10.1038/ng.3635
    [76] Boque-Sastre R, Soler M, Oliveira-Mateos C, et al. (2015) Head-to-head antisense transcription and R-loop formation promotes transcriptional activation. P Natl Acad Sci USA 112: 5785–5790. doi: 10.1073/pnas.1421197112
  • This article has been cited by:

    1. Shuqi Zhai, Qinglong Wang, Ting Yu, Fuzzy optimal harvesting of a prey-predator model in the presence of toxicity with prey refuge under imprecise parameters, 2022, 19, 1551-0018, 11983, 10.3934/mbe.2022558
    2. Yuan Tian, Chunxue Li, Jing Liu, Complex dynamics and optimal harvesting strategy of competitive harvesting models with interval-valued imprecise parameters, 2023, 167, 09600779, 113084, 10.1016/j.chaos.2022.113084
    3. I. Sukarsih, A. K. Supriatna, E. Carnia, N. Anggriani, A Runge-Kutta numerical scheme applied in solving predator-prey fuzzy model with Holling type II functional response, 2023, 9, 2297-4687, 10.3389/fams.2023.1096167
    4. Xuyang Cao, Qinglong Wang, Jie Liu, Hopf bifurcation in a predator-prey model under fuzzy parameters involving prey refuge and fear effects, 2024, 9, 2473-6988, 23945, 10.3934/math.20241164
    5. Christian Cortés García, Population dynamics in a Leslie–Gower predator–prey model with predator harvesting at high densities, 2024, 0170-4214, 10.1002/mma.10359
    6. Yuan Tian, Hua Guo, Wenyu Shen, Xinrui Yan, Jie Zheng, Kaibiao Sun, Dynamic analysis and validation of a prey-predator model based on fish harvesting and discontinuous prey refuge effect in uncertain environments, 2025, 33, 2688-1594, 973, 10.3934/era.2025044
  • Reader Comments
  • © 2018 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(5084) PDF downloads(952) Cited by(1)

Figures and Tables

Figures(1)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog