Research article Special Issues

Meta-analysis identifying epithelial-derived transcriptomes predicts poor clinical outcome and immune infiltrations in ovarian cancer


  • Received: 21 June 2021 Accepted: 19 July 2021 Published: 30 July 2021
  • Background 

    Previous studies revealed that the epithelial component is associated with the modulation of the ovarian tumor microenvironment (TME). However, the identification of key transcriptional signatures of laser capture microdissected human ovarian cancer epithelia remains lacking.

    Methods 

    We identified the differentially expressed transcriptional signatures of human ovarian cancer epithelia by meta-analysis of GSE14407, GSE2765, GSE38666, GSE40595, and GSE54388. Then we investigated the enrichment of KEGG pathways that are associated with epithelia-derived transcriptomes. Finally, we investigated the correlation of key epithelia-hub genes with the survival prognosis and immune infiltrations. Finally, we investigated the genetic alterations of key prognostic hub genes and their diagnostic efficacy in ovarian cancer epithelia.

    Results 

    We identified 1339 differentially expressed genes (DEGs) in ovarian cancer epithelia including 541upregulated and 798 downregulated genes. We identified 21 (such as E2F4, FOXM1, TFDP1, E2F1, and SIN3A) and 11 (such as JUN, DDX4, FOSL1, NOC2L, and HMGA1) master transcriptional regulators (MTRs) that are interacted with upregulated and the downregulated genes in ovarian tumor epithelium, respectively. The STRING-based analysis identified hub genes (such as CDK1, CCNB1, AURKA, CDC20, and CCNA2) in ovarian cancer epithelia. The significant clusters of identified hub genes are associated with the enrichment of KEGG pathways including cell cycle, DNA replication, cytokine-cytokine receptor interaction, pathways in cancer, and focal adhesion. The upregulation of SCNN1A and CDCA3 and the downregulation of SOX6 are correlated with a shorter survival prognosis in ovarian cancer (OV). The expression level of SOX6 is negatively correlated with immune score and positively correlated with tumor purity in OV. Moreover, SOX6 is negatively correlated with the infiltration of TILs, CD8+ T cells, CD4+ Regulatory T cells, cytolytic activity, T cell activation, pDC, neutrophils, and macrophages in OV. Also, SOX6 is negatively correlated with various immune markers including CD8A, PRF1, GZMA, GZMB, NKG7, CCL3, and CCL4, indicating the immune regulatory efficiency of SOX6 in the TME of OV. Furthermore, SCNN1A, CDCA3, and SOX6 genes are genetically altered in OV and the expression levels of SCNN1A and SOX6 genes showed diagnostic efficacy in ovarian cancer epithelia.

    Conclusions 

    The identified ovarian cancer epithelial-derived key transcriptional signatures are significantly correlated with survival prognosis and immune infiltrations, and may provide new insight into the diagnosis and treatment of epithelial ovarian cancer.

    Citation: Dong-feng Li, Aisikeer Tulahong, Md. Nazim Uddin, Huan Zhao, Hua Zhang. Meta-analysis identifying epithelial-derived transcriptomes predicts poor clinical outcome and immune infiltrations in ovarian cancer[J]. Mathematical Biosciences and Engineering, 2021, 18(5): 6527-6551. doi: 10.3934/mbe.2021324

    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
  • Background 

    Previous studies revealed that the epithelial component is associated with the modulation of the ovarian tumor microenvironment (TME). However, the identification of key transcriptional signatures of laser capture microdissected human ovarian cancer epithelia remains lacking.

    Methods 

    We identified the differentially expressed transcriptional signatures of human ovarian cancer epithelia by meta-analysis of GSE14407, GSE2765, GSE38666, GSE40595, and GSE54388. Then we investigated the enrichment of KEGG pathways that are associated with epithelia-derived transcriptomes. Finally, we investigated the correlation of key epithelia-hub genes with the survival prognosis and immune infiltrations. Finally, we investigated the genetic alterations of key prognostic hub genes and their diagnostic efficacy in ovarian cancer epithelia.

    Results 

    We identified 1339 differentially expressed genes (DEGs) in ovarian cancer epithelia including 541upregulated and 798 downregulated genes. We identified 21 (such as E2F4, FOXM1, TFDP1, E2F1, and SIN3A) and 11 (such as JUN, DDX4, FOSL1, NOC2L, and HMGA1) master transcriptional regulators (MTRs) that are interacted with upregulated and the downregulated genes in ovarian tumor epithelium, respectively. The STRING-based analysis identified hub genes (such as CDK1, CCNB1, AURKA, CDC20, and CCNA2) in ovarian cancer epithelia. The significant clusters of identified hub genes are associated with the enrichment of KEGG pathways including cell cycle, DNA replication, cytokine-cytokine receptor interaction, pathways in cancer, and focal adhesion. The upregulation of SCNN1A and CDCA3 and the downregulation of SOX6 are correlated with a shorter survival prognosis in ovarian cancer (OV). The expression level of SOX6 is negatively correlated with immune score and positively correlated with tumor purity in OV. Moreover, SOX6 is negatively correlated with the infiltration of TILs, CD8+ T cells, CD4+ Regulatory T cells, cytolytic activity, T cell activation, pDC, neutrophils, and macrophages in OV. Also, SOX6 is negatively correlated with various immune markers including CD8A, PRF1, GZMA, GZMB, NKG7, CCL3, and CCL4, indicating the immune regulatory efficiency of SOX6 in the TME of OV. Furthermore, SCNN1A, CDCA3, and SOX6 genes are genetically altered in OV and the expression levels of SCNN1A and SOX6 genes showed diagnostic efficacy in ovarian cancer epithelia.

    Conclusions 

    The identified ovarian cancer epithelial-derived key transcriptional signatures are significantly correlated with survival prognosis and immune infiltrations, and may provide new insight into the diagnosis and treatment of epithelial ovarian cancer.



    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 ˜A in a universe of discourse X is defined as the following set of pairs ˜A={(x,μ˜A(x)):xX}. The mapping μ˜A:X[0,1] is called the membership function of the fuzzy set ˜A and μ˜A is called the membership value or degree of membership of xX in the fuzzy set ˜A.

    Definition 2.2. [54] α-cut of fuzzy set: The α-cut of a fuzzy set ˜A is a crisp set and it is defined by Aα={x:μ˜A(x)α},α(0,1]. For α=0 the support of ˜A is defined as A0=Supp(˜A)=¯{xR,μ˜A(x)>0}.

    Definition 2.3. [32] Convex fuzzy set: A convex fuzzy set ˜A is a fuzzy set on a continuous universe such that for all α, Aα 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) ˜A(a1,a2,a3) is fuzzy set of the real line R characterized by the membership function μ˜A:R[0,1] as follows

    μ˜A={xa1a2a1  if  a1xa2,a3xa3a2  if  a2xa3,0  otherwise.

    Therefore, the α-cut of triangular fuzzy number is a bounded and closed internal [AL(α),AR(α)], where AL(α)=inf{x:μ˜A(x)α}=a1+α(a2a1) and AR(α)=sup{x:μ˜A(x)α}=a3α(a3a2).

    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 ˜AE, then the α-level set Aα is a bounded closed interval which can be expressed as Aα=[AαL,AαR], where AαL=inf{x:μ˜A(x)α} and AαR=sup{x:μ˜A(x)α}.

    Lemma 2.1. [57] If ˜A,˜BE, then for α(0,1],

    (1) [˜A+˜B]α=[AαL+BαL,AαR+BαR],

    (2) [˜A˜B]α=[min{AαLBαL,AαLBαR,AαRBαL,AαRBαR},max{AαLBαL,AαLBαR,AαRBαL,AαRBαR}],

    (3) [˜A˜B]α=[AαLBαR,AαRBαL].

    Lemma 2.2. [58] [AαL,AαR],0<α1 is a family of nonempty intervals. If

    (1) [AαL,AαR][AβL,AβR] for 0<αβ and

    (2) [limkAαkL,limkAαkR]=[AαL,AαR],

    whenever non-decreasing sequence αk converges to α(0,1]; then, the family [AαL,AαR], 0<α1, represents the α-level sets of a fuzzy number ˜A in E. Conversely, if [AαL,AαR], 0<α1 are the α-level sets of a fuzzy number ˜AE, then the conditions (1) and (2) hold true.

    In weighted sum method [59], a simple utility function is in the form of wigi for ith objective, where wi stands for the weight of ith objective. And the total utility function U is expressed as follows

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

    where wi>0 and liwi=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 x1(0)>0 and x2(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, q1,q2 and E1,E2 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 ˜r,˜c,˜d,˜s and ˜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 ˜r while predator decreases with fuzzy mortality ˜d for lack of prey. ˜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 ˜c(1m)x, and ˜e represents the fuzzy growth rate of predator attacking on prey.

    Applying α-level to cut triangular fuzzy numbers ˜r,˜c,˜d,˜s and ˜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 w1 and w2 express as weight satisfying w1+w2=1, and w1,w20. 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 (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,ξ), where 0<ξ<+. 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 lim suptx(t)A1KA2k1. From the second differential equation of system (3.5), one has

    dydtk1B3(1m)yB2y2,

    which can be written as

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

    It is easy to see lim supty(t)k1B3(1m)B2. 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: S0=(0,0).

    (2) Axial equilibrium: S1=(K(A1q1E1)A2,0) exists if A1>q1E1.

    (3) Axial equilibrium: S2=(0,q2E2+B1B2) does not exist since q2E2+B1B2<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 B3(1m)x>B1+q2E2 and A1B2+A3(1m)(B1+q2E2)>B2q1E1 are satisfied.

    On account of the following theorem, the local stability of equilibria S0,S1,S for system (3.5) is discussed.

    Theorem 5.1. The following conclusions hold:

    (1) Trivial equilibrium S0 becomes stable node when A1<q1E1 or saddle point when A1>q1E1.

    (2) Axial equilibrium S1 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)

    (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 M0 at S0 is represented as

    M0=(A1q1E100B1q2E2).

    According to the characteristic equation about the above matrix

    (λA1+q1E1)(λ+B1+q2E2)=0,

    one has λ1=A1q1E1,  λ2=B1q2E2, hence A0 is stable node when A1<q1E1 or saddle point when A1>q1E1.

    (2) The Jacobian matrix M1 at S1 is depicted as

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

    In accordance with the characteristic equation of the above matrix

    (λ+A1q1E1)(λ+B1+q2E2KB3(1m)(A1q1E1)A2)=0,

    we have λ1=(A1q1E1) and λ2=(B1+q2E2)+KB3(1m)(A1q1E1)A2. Obviously, S1 exists if A1>q1E1 holds, hence S1 is a stable node when A1>q1E1 and A2(B1+q2E2)>KB3(1m)(A1q1E1) or saddle point when A1>q1E1 and A2(B1+q2E2)<KB3(1m)(A1q1E1).

    (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., λ1+λ2=(B2y+A2xK)<0 and the product of the roots of Eq (5.6) is positive, i.e., λ1λ2=A3B3(1m)2xy>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 B3(1m)x>B1+q2E2 and A1B2+A3(1m)(B1+q2E2)>B2q1E1 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=A3B3, it yields to

    dVdt=A2K(xx)2A3B2B3(yy)2<0.

    Since dVdt 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 dxdt0 and dydt0. 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, c1,c2 denote unit capture cost for prey x and predator y, respectively. And p1,p2 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 N1=(p1q1xc1)E1,N2=(p2q2yc2)E2. N1 and N2 separately represent the profit function of prey x and predator y. And the bionomic equilibrium (x,y,E1,E2) 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: E1=0. One has

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

    Clearly, (x,y,E1,E2)=(0,c2p2q2,0,p2q2B1+B2c2p2q22) satisfies Eq (6.3), we have to omit the bionomic equilibrium owing to E2<0.

    Then solving the second equation for a nonzero x yields

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

    which exists if satisfying A1p2q2>A3(1m)c2. 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,y,0,E2).

    Case Ⅱ Boundary bionomic equilibria with no catching species y: E2=0. It follows from

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

    that (x,y,E1,E2)=(c1p1q1,0,A1p1q1KA2c1p1q21K,0) and we get a nonegative semi-trivial bionomic equilibrium (x,0,E1,0) with condition A1p1q1K>A2c1.

    Then solving the third equation for a nonzero y yields

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

    satisfying B3c1(1m)>B1p1q1. 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,y,E1,0).

    Case Ⅲ Boundary bionomic equilibria with no catching species x and y: E1=0 and E2=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 E1=E2=0.

    Case Ⅳ Positive bionomic equilibrium with effort fishing on both species x and y: E1>0 and E2>0.

    Solving the third equation of Eq (6.2), we obtain x=c1p1q1,y=c2p2q2. 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,y,E1,E2).

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

    Theorem 6.1. The following conclusions hold:

    (1) The boundary bionomic equilibrium (x,y,0,E2) exists if A1p2q2>A3(1m)c2 and Eq (6.6) are satisfied, where y=c2p2q2 and x and E2 are defined in Eqs (6.4) and (6.5), respectively.

    (2) The boundary bionomic equilibrium (x,0,E1,0) exists if A1p1q1K>A2c1 holds, where x=c1p1q1 and E1=A1p1q1KA2c1p1q21K.

    (3) The boundary bionomic equilibrium (x,y,E1,0) exists if B3c1(1m)>B1p1q1 and Eq (6.10) are satisfied, where x=c1p1q1 and y and E1 are defined in Eqs (6.8) and (6.9), respectively.

    (4) The boundary bionomic equilibrium (x,y,0,0) is equal to the biological equilibrium with E1=E2=0 described in Section 5.

    (5) The positive bionomic equilibrium (x,y,E1,E2) exists if Eq (6.12) holds, where x=c1p1q1, y=c2p2q2 and E1 and E2 are defined in Eq (6.11).

    This section studies the optimal harvesting strategy with fuzzy net discount rate of inflation. Previously, denote ˜k and ˜r the inflation and discount rates and they are considered as fuzzy number essentially. The net discount rate of inflation ˜δ is the difference of the above two fuzzy numbers, that is ˜δ=˜r˜k, and it can be regarded as triangular fuzzy number, that is ˜δ=(δ1,δ2,δ3). Our target is to maximize the value ˜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 Ei(t) subject to 0Ei(t)Emaxi,i=1,2. According to Maity and Maiti [61], Sadhukhan et al. [52] and Pal and Mahapatra [23], using α to cut triangular fuzzy number ˜δ=(δ1,δ2,δ3), that is it regarded as an interval number [δL,δR], the optimization problem Maximize ˜J(E1,E2) is equal to Maximize [JL(E1,E2),JR(E1,E2)], 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 w1,w20 are two weights such that w1+w2=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 λi=λ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 λ1 in Eq (7.7) to get a reduced differential equation for λ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 Ci(i=1,2) are arbitrary constants and mi(i=1,2) are the roots of the auxiliary equations a0m2+a1m+a2=0 and

    NL=a0δ2Ra1δR+a20,   NR=a0δ2La1δL+a20.

    From Eq (7.9), λ2 is bounded if and only if mi<0(i=1,2) or Ci=0. Since it is difficult to check whether mi<0, we consider Ci=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 λieδLt(i=1,2) of the two species satisfy the transversality condition at [18], i.e., they are bounded as t. The Hamiltonian in Eq (7.4) must be maximized for Ei[0,Emaxi]. Suppose that the optimal equilibrium does not appear at either Ei=0 or Ei=Emaxi, 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 λ1 and λ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(δRδL)t in Eq (7.17), we consider the following two cases:

    Case Ⅰ If PLc1w1 and QLc2w1, then Eq (7.17) is equal to

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

    Case Ⅱ If PL=c1w1 or QL=c2w1, we have

    PLc1w1=PRc1w2=0,  or  QLc2w1=QRc2w2=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 λ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 λ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 E1 and E2 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 N1=(p1q1xc1)E1>0 and N2=(p2q2yc2)E2>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 E1 and E2 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˜δ,y=y˜δ and the optimal harvesting efforts E1=E1˜δ,E2=E2˜δ.

    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: ˜r=(2.00,2.20,2.40), ˜c=(1.20,1.35,1.50), ˜d=(0.30,0.40,0.50), ˜s=(0.06,0.07,0.08), ˜e=(0.60,0.70,0.80), K=5 and m=0.1.

    We discuss equilibria S0, S1 and S of the predator-prey model under different values of α,w1 and w2, which are showed from Table 1 to Table 3.

    Table 1.  The trivial equilibrium S0 for q1=0.2,q2=0.2,E1=15 and E2=10.
    w1 w2 S0 at α=0 S0 at α=0.3 S0 at α=0.6 S0 at α=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 S1 for q1=0.2,q2=0.2,E1=5 and E2=15.
    w1 w2 S1 at α=0 S1 at α=0.3 S1 at α=0.6 S1 at α=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 q1=0.2,q2=0.2,E1=5 and E2=2.
    w1 w2 S at α=0 S at α=0.3 S at α=0.6 S at α=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 w1, w2 and α, the values of prey x and predator y stay at zero; for the axial equilibrium, the value of prey x decreases under the same α with increasing w1, predator y is invariant in zero; for the interior equilibrium, the value of prey x increases under the same α with increasing w1, and predator y decreases under the same α with increasing w1.

    Figures 13 are time series of prey x and predator y with initial values (2,6) for different w1, w2 and α, 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 q1=0.2,q2=0.2,E1=15,E2=10 and different w1,w2 and α, respectively, t[0,100].
    Figure 2.  (a)–(f) Time series diagram of prey x and predator y with initial values (2,6) for q1=0.2,q2=0.2,E1=5,E2=15 and different w1,w2 and α, respectively, t[0,100].
    Figure 3.  (a)–(f) Time series diagram of prey x and predator y with initial values (2,6) for q1=0.2,q2=0.2,E1=5,E2=2 and different w1,w2 and α, respectively, t[0,100].

    The phase portrait of prey x and predator y corresponding to stabilities of different equilibria for all kinds of combinations of w1, w2 and α 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 q1=0.2,q2=0.2,E1=15,E2=10 and different w1, w2 and α, respectively, t[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 q1=0.2,q2=0.2,E1=5,E2=15 and different w1, w2 and α, respectively, t[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 q1=0.2,q2=0.2,E1=5,E2=2 and different w1, w2 and α, respectively, t[0,100].

    Dynamical behavior of prey x and predator y under interior steady state level with respect to α for initial condition (x(0),y(0))=(2.5,6), and q1=0.2,q2=0.2,E1=5,E2=2 are presented in Figure 7.

    Figure 7.  (a)–(f) Dynamical behavior of prey x and predator y with respect to α for initial condition (x(0),y(0))=(2.5,6), and q1=0.2,q2=0.2,E1=5,E2=2.

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

    In the mean time, the dynamical behavior of prey x and predator y about w1 and w2 under initial condition (x(0),y(0))=(2.5,6) for α=0,0.3,0.6 and α=0.9, and q1=0.2,q2=0.2,E1=5,E2=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 α with increasing w1, and predator y decreases under fixed α with increasing w1.

    Figure 8.  (a)–(d) Dynamical behavior of prey x and predator y with respect to w1 and w2 under initial condition (x(0),y(0))=(2.5,6) for α=0,0.3,0.6 and α=0.9 respectively, and q1=0.2, q2=0.2, E1=5, E2=2.

    Example 2 Consider the parameter values as follows: ˜r=(1.50,1.55,1.60), ˜c=(0.250,0.275,0.300), ˜d=(0.450,0.475,0.500), ˜s=(0.150,0.175,0.200), ˜e=(1.300,1.325,1.350), K=10,m=0.1,q1=0.92,q2=0.95,p1=20,p2=25,c1=30 and c2=15.

    The bionomic equilibria for different combinations of w1, w2 and α are emerged in Table 4. It is observed that x and y are invariable with respect to w1, w2 and α, and E1 and E2 are decreasing with fixed α and increasing w1. As well as, together with the value of α increasing, E1 and E2 decreases for w1=0,w2=1 and w1=0.2,w2=0.8 and w1=0.4,w2=0.6. Nevertheless, in the cases w1=0.6,w2=0.4 and w1=0.8,w2=0.2 and w1=1,w2=0, as the value of α increases, E1 and E2 are increasing.

    Table 4.  Positive bionomic equilibrium for different w1,w2 and α.
    w1 w2 α=0 α=0.3 α=0.6 α=0.9
    (x,y) (E1,E2) (x,y) (E1,E2) (x,y) (E1,E2) (x,y) (E1,E2)
    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: ˜r=(2.80,2.83,2.85),˜c=(2.45,2.55,2.65),˜d=(0.012,0.013,0.015),˜s=(0.008,0.009,0.010),˜e=(0.20,0.21,0.22),K=10,m=0.1,q1=0.45,q2=0.55,p1=30,p2=25,c1=25,c2=12 and ˜δ=(0.8,0.9,1.0).

    The optimal equilibrium and harvesting effort for different combinations w1,w2,α and w1,w2,α are appeared in Tables 59, respectively. For different combinations w1,w2,α and w1,w2,α, we see the existence of optimal equilibrium.

    Table 5.  Optimal equilibrium and optimal harvesting effort for w1=0,w2=1,α=0.
    w1 w2 α=0 α=0.3 α=0.6 α=0.9
    (x˜δ,y˜δ) (E1˜δ,E2˜δ) (x˜δ,y˜δ) (E1˜δ,E2˜δ) (x˜δ,y˜δ) (E1˜δ,E2˜δ) (x˜δ,y˜δ) (E1˜δ,E2˜δ)
    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 w1=0,w2=1,α=0.5.
    w1 w2 α=0 α=0.3 α=0.6 α=0.9
    (x˜δ,y˜δ) (E1˜δ,E2˜δ) (x˜δ,y˜δ) (E1˜δ,E2˜δ) (x˜δ,y˜δ) (E1˜δ,E2˜δ) (x˜δ,y˜δ) (E1˜δ,E2˜δ)
    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 w1=0,w2=1,α=0.9.
    w1 w2 α=0 α=0.3 α=0.6 α=0.9
    (x˜δ,y˜δ) (E1˜δ,E2˜δ) (x˜δ,y˜δ) (E1˜δ,E2˜δ) (x˜δ,y˜δ) (E1˜δ,E2˜δ) (x˜δ,y˜δ) (E1˜δ,E2˜δ)
    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 w1=0.5,w2=0.5,α=0.9.
    w1 w2 α=0 α=0.3 α=0.6 α=0.9
    (x˜δ,y˜δ) (E1˜δ,E2˜δ) (x˜δ,y˜δ) (E1˜δ,E2˜δ) (x˜δ,y˜δ) (E1˜δ,E2˜δ) (x˜δ,y˜δ) (E1˜δ,E2˜δ)
    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 w1=1,w2=0,α=0.9.
    w1 w2 α=0 α=0.3 α=0.6 α=0.9
    (x˜δ,y˜δ) (E1˜δ,E2˜δ) (x˜δ,y˜δ) (E1˜δ,E2˜δ) (x˜δ,y˜δ) (E1˜δ,E2˜δ) (x˜δ,y˜δ) (E1˜δ,E2˜δ)
    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 α=0 it is observed that prey x˜δ and predator y˜δ are invariant, optimal harvesting effort of the prey E1˜δ decreases and then increases, but the optimal harvesting effort of the predator E2˜δ increases and then decreases with increasing w1; when α=0.3, x˜δ decreases and then increases, y˜δ is invariant, E1˜δ decreases and then increases, and E2˜δ increases and then decreases with increasing w1; when α=0.6, x˜δ increases and then decreases, y˜δ is invariant, E1˜δ fluctuates between 0.0907 and 0.0908, and E2˜δ is invariant with increasing w1; when α=0.9, x˜δ and y˜δ are invariant, E1˜δ decreases and then increases, and E2˜δ is invariant with increasing w1. Besides, when fixed w1 and w2, x˜δ increases, y˜δ decreases, E1˜δ increases, and E2˜δ decreases with the increase of α.

    In Table 6, when fixed α, x˜δ, y˜δ, E1˜δ and E2˜δ are invariant with the increase of w1; when fixed w1 and w2, x˜δ increases, y˜δ decreases, E1˜δ increases, and E2˜δ increases with the increase of α.

    In Table 7, when fixed α, x˜δ, y˜δ, E1˜δ and E2˜δ are invariant with the increase of w1; when fixed w1 and w2, x˜δ increases, y˜δ is invariant, E1˜δ decreases, and E2˜δ increases with the increase of α.

    In Table 8, when fixed α, x˜δ, y˜δ, E1˜δ and E2˜δ are invariant with the increase of w1; when fixed w1 and w2, x˜δ increases, y˜δ increases, E1˜δ decreases, and E2˜δ increases with the increase of α.

    In Table 9, when α=0, x˜δ increases and then decreases and y˜δ is invariant, E1˜δ increases and then decreases, but E2˜δ decreases and then increases with increasing w1; when α=0.3 or α=0.6, we can see that x˜δ increases and then decreases, y˜δ, E1˜δ and E2˜δ are invariant with increasing w1; when α=0.9, x˜δ, y˜δ and E1˜δ are invariant, E2˜δ decreases and then increases with increasing w1. When fixed w1 and w2, x˜δ decreases, y˜δ is invariant, E1˜δ is almost invariant, and E2˜δ increases with the increase of α.

    From Tables 57, it is easy to see that when fixed w1,w2 and w1,w2,α, 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 α.

    Based on Tables 79, we observed that when fixed α and w1,w2,α, 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 w1.

    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 w1,w2 and α. Dynamical behaviors of prey and predator about α or w1 and w2, 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 ˜J. The optimal steady state solution is displayed in Tables 59. We obtain optimal equilibrium and optimal harvesting effort as respect to w1,w2 and α under five cases including w1=0, w2=1, α=0 and w1=0, w2=1, α=0.5 and w1=0, w2=1, α=0.9 and w1=0.5, w2=0.5, α=0.9 and w1=1, w2=0, α=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 α-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] L. M. Hurwitz, P. F. Pinsky, B. Trabert, General population screening for ovarian cancer, Lancet, 397 (2021), 2128-2130. doi: 10.1016/S0140-6736(21)01061-8
    [2] R. L. Siegel, K. D. Miller, H. E. Fuchs, A. Jemal, Cancer statistics, 2021, CA Cancer J. Clin., 71 (2021), 7-33. doi: 10.3322/caac.21654
    [3] N. J. Bowen, L. D. Walker, L. V. Matyunina, S. Logani, K. A. Totten, B. B. Benigno, et al., Gene expression profiling supports the hypothesis that human ovarian surface epithelia are multipotent and capable of serving as ovarian cancer initiating cells, BMC Med. Genomics., 2 (2009), 71. doi: 10.1186/1755-8794-2-71
    [4] T. L. Yeung, C. S. Leung, K. K. Wong, G. Samimi, M. S. Thompson, J. Liu, et al., TGF-β modulates ovarian cancer invasion by upregulating CAF-derived versican in the tumor microenvironment, Cancer Res., 73 (2013), 5016-5028. doi: 10.1158/0008-5472.CAN-13-0023
    [5] J. Farley, L. L. Ozbun, M. J. Birrer, Genomic analysis of epithelial ovarian cancer, Cell Res., 18 (2008), 538-548. doi: 10.1038/cr.2008.52
    [6] M. Ezzati, A. Abdullah, A. Shariftabrizi, J. Hou, M. Kopf, J. K. Stedman, et al., Recent advancements in prognostic factors of epithelial ovarian carcinoma, Int. Sch. Res. Not., 2014 (2014), e953509.
    [7] M. L. Drakes, P. J. Stiff, Regulation of ovarian cancer prognosis by immune cells in the tumor microenvironment, Cancers, 10 (2018), 302. doi: 10.3390/cancers10090302
    [8] T. Gui, C. Yao, B. Jia, K. Shen, Identification and analysis of genes associated with epithelial ovarian cancer by integrated bioinformatics methods, PloS One, 16 (2021), e0253136. doi: 10.1371/journal.pone.0253136
    [9] J. Liu, H. Meng, S. Li, Y. Shen, H. Wang, W. Shan, et al., Identification of potential biomarkers in association with progression and prognosis in epithelial ovarian cancer by integrated bioinformatics analysis, Front. Genet., 10 (2019), 1031. doi: 10.3389/fgene.2019.01031
    [10] D. Matei, T. G. Graeber, R. L. Baldwin, B. Y. Karlan, J. Rao, D. D. Chang, Gene expression in epithelial ovarian carcinoma, Oncogene, 21 (2002), 6289-6298. doi: 10.1038/sj.onc.1205785
    [11] P. Israelsson, E. Dehlin, I. Nagaev, E. Lundin, U. Ottander, L. Mincheva-Nilsson, Cytokine mRNA and protein expression by cell cultures of epithelial ovarian cancer-methodological considerations on the choice of analytical method for cytokine analyses, Am. J. Reprod. Immunol., 84 (2020), e13249.
    [12] K. P. Prahm, C. K. Hø gdall, M. A. Karlsen, I. J. Christensen, G. W. Novotny, E. Hø gdall, MicroRNA characteristics in epithelial ovarian cancer, PloS One, 16 (2021), e0252401. doi: 10.1371/journal.pone.0252401
    [13] A. A. Alshamrani, Roles of microRNAs in ovarian cancer tumorigenesis: two decades later, what have we learned?, Front. Oncol., 10 (2020), 1084.
    [14] S. Udhaya Kumar, D. Thirumal Kumar, R. Bithia, S. Sankar, R. Magesh, M. Sidenna, et al., Analysis of differentially expressed genes and molecular pathways in familial hypercholesterolemia involved in atherosclerosis: a systematic and bioinformatics approach, Front. Genet., 11 (2020), 734. doi: 10.3389/fgene.2020.00734
    [15] J. Wan, S. Jiang, Y. Jiang, W. Ma, X. Wang, Z. He, et al., Corrigendum: data mining and expression analysis of differential lncRNA ADAMTS9-AS1 in prostate cancer, Front. Genet., 11 (2020), 361. doi: 10.3389/fgene.2020.00361
    [16] U. K. S., B. Rajan, T. K. D., A. P. V., T. Abunada, S. Younes, et al., Involvement of essential signaling cascades and analysis of gene networks in diabesity, Genes, 11 (2020), 1256. doi: 10.3390/genes11111256
    [17] D. Fu, B. Zhang, L. Yang, S. Huang, W. Xin, Development of an immune-related risk signature for predicting prognosis in lung squamous cell carcinoma, Front. Genet., 11 (2020), 978. doi: 10.3389/fgene.2020.00978
    [18] E. R. King, C. S. Tung, Y. T. M. Tsang, Z. Zu, G. T. M. Lok, M. T. Deavers, et al, The anterior gradient homolog 3 (AGR3) gene is associated with differentiation and survival in ovarian cancer, Am. J. Surg. Pathol., 35 (2011), 904-912.
    [19] C. Huang, E. A. Clayton, L. V. Matyunina, L. D. McDonald, B. B. Benigno, F. Vannberg, et al., Machine learning predicts individual cancer patient responses to therapeutic drugs with high accuracy, Sci. Rep., 8 (2018), 16444. doi: 10.1038/s41598-018-34753-5
    [20] L. N. Lili, L. V. Matyunina, L. D. Walker, B. B. Benigno, J. F. McDonald, Molecular profiling predicts the existence of two functionally distinct classes of ovarian cancer stroma, BioMed Res. Int., 2013 (2013), 846387.
    [21] T. L. Yeung, C. S. Leung, K. K. Wong, A. Gutierrez-Hartmann, J. Kwong, D. M. Gershenson, et al., ELF3 is a negative regulator of epithelial-mesenchymal transition in ovarian cancer cells, Oncotarget, 8 (2017), 16951-16963. doi: 10.18632/oncotarget.15208
    [22] J. Xia, E. E. Gill, R. E. W. Hancock, NetworkAnalyst for statistical, visual and network-based meta-analysis of gene expression data, Nat. Protoc., 10 (2015), 823-844.
    [23] W. E. Johnson, C. Li, A. Rabinovic, Adjusting batch effects in microarray expression data using empirical Bayes methods, Biostat. Oxf. Engl., 8 (2007), 118-127.
    [24] W. G. Cochran, The combination of estimates from different experiments, Biometrics, 10 (1954), 101-129. doi: 10.2307/3001666
    [25] Y. Benjamini, Y. Hochberg, Controlling the false discovery rate: a practical and powerful approach to multiple testing, J. R. Stat. Soc. Ser. B Methodol., 57 (1995), 289-300.
    [26] P. Shannon, A. Markiel, O. Ozier, N. S. Baliga, J. T. Wang, D. Ramage, et al., Cytoscape: a software environment for integrated models of biomolecular interaction networks, Genome Res., 13 (2003), 2498-2504. doi: 10.1101/gr.1239303
    [27] R. Janky, A. Verfaillie, H. Imrichová, B. Van de Sande, L. Standaert, V. Christiaens, et al., iRegulon: from a gene list to a gene regulatory network using large motif and track collections, PLoS Comput. Biol., 10 (2014), e1003731.
    [28] J. Wang, Y. Wu, M. N. Uddin, R. Chen, J. P. Hao, Identification of potential key genes and regulatory markers in essential thrombocythemia through integrated bioinformatics analysis and clinical validation, Pharmacogenomics Pers. Med., 14 (2021), 767-784.
    [29] A. Subramanian, P. Tamayo, V. K. Mootha, S. Mukherjee, B. L. Ebert, M. A. Gillette, et al., Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles, Proc. Natl. Acad. Sci., 102 (2005), 15545-15550. doi: 10.1073/pnas.0506580102
    [30] M. Kanehisa, M. Furumichi, M. Tanabe, Y. Sato, K. Morishima, KEGG: new perspectives on genomes, pathways, diseases and drugs, Nucleic Acids Res., 45 (2017), D353-D361.
    [31] D. Szklarczyk, A. L. Gable, D. Lyon, A. Junge, S. Wyder, J. Huerta-Cepas, et al., STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets, Nucleic Acids Res., 47 (2019), D607-D613.
    [32] C. H. Chin, S. H. Chen, H. H. Wu, C. W. Ho, M. T. Ko, C. Y. Lin, cytoHubba: identifying hub objects and sub-networks from complex interactome, BMC Syst. Biol., 8 (2014), S11.
    [33] G. D. Bader, C. W. Hogue, An automated method for finding molecular complexes in large protein interaction networks, BMC Bioinformatics, 4 (2003), 2. doi: 10.1186/1471-2105-4-2
    [34] T. Therneau, A package for survival analysis in R, 2021. Available from: https://cran.r-project.org/web/packages/survival/vignettes/survival.pdf.
    [35] K. Yoshihara, M. Shahmoradgoli, E. Martínez, R. Vegesna, H. Kim, W. Torres-Garcia, et al., Inferring tumour purity and stromal and immune cell admixture from expression data, Nat. Commun., 4 (2013), 2612. doi: 10.1038/ncomms3612
    [36] S. Hä nzelmann, R. Castelo, J. Guinney, GSVA: gene set variation analysis for microarray and RNA-seq data, BMC Bioinformatics, 14 (2013), 7. doi: 10.1186/1471-2105-14-7
    [37] I. Tirosh, B. Izar, S. M. Prakadan, M. H. Wadsworth, D. Treacy, J. J. Trombetta, et al., Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq, Science, 352 (2016), 189-196. doi: 10.1126/science.aad0501
    [38] D. Aran, Z. Hu, A. J. Butte, xCell: digitally portraying the tissue cellular heterogeneity landscape, Genome Biol., 18 (2017), 220. doi: 10.1186/s13059-017-1349-1
    [39] T. Davoli, H. Uno, E. C. Wooten, S. J. Elledge, Tumor aneuploidy correlates with markers of immune evasion and with reduced response to immunotherapy, Science, 355 (2017), eaaf8399.
    [40] Z. Liu, M. Li, Z. Jiang, X. Wang, A comprehensive immunologic portrait of triple-negative breast cancer, Transl. Oncol., 11 (2018), 311-329. doi: 10.1016/j.tranon.2018.01.011
    [41] E. Cerami1, J. Gao, U. Dogrusoz, B. E. Gross, S. O. Sumer, B. A. Aksoy, et al., The cBio cancer genomics portal: an open platform for exploring multidimensional cancer genomics data, Cancer Discov., 2 (2012), 401-404. doi: 10.1158/2159-8290.CD-12-0095
    [42] X. Robin, N. Turck, A. Hainard, N. Tiberti, F. Lisacek, J. C. Sanchez, et al., pROC: an open-source package for R and S+ to analyze and compare ROC curves, BMC Bioinformatics, 12 (2011), 77. doi: 10.1186/1471-2105-12-77
    [43] Q. Yang, R. Wang, B. Wei, C. Peng, L. Wang, G. Hu, et al., Candidate biomarkers and molecular mechanism investigation for glioblastoma multiforme utilizing WGCNA, BioMed Res. Int., 2018 (2018), e4246703.
    [44] M. N. Uddin, R. Akter, M. Li, Z. Abdelrahman, Expression of SARS-COV-2 cell receptor gene ACE2 is associated with immunosuppression and metabolic reprogramming in lung adenocarcinoma based on bioinformatics analyses of gene expression profiles, Chem. Biol. Interact., 335 (2021), 109370. doi: 10.1016/j.cbi.2021.109370
    [45] V. Gómez-Rubio, ggplot2 - elegant graphics for data analysis, J. Stat. Softw., 77 (2017), 1-3.
    [46] L. Zhao, Y. Li, Z. Zhang, J. Zou, J. Li, R. Wei, et al., Meta-analysis based gene expression profiling reveals functional genes in ovarian cancer, Biosci. Rep., 40 (2020), BSR20202911.
    [47] D. Chaves-Moreira, M. A. Mitchell, C. Arruza, P. Rawat, S. Sidoli, R. Nameki, et al., PAX8 orchestrates an angiogenic program through interaction with SOX17, preprint, BioRxiv, 2020.09.09.290387.
    [48] X. Liu, Y. Gao, Y. Lu, J. Zhang, L. Li, F. Yin, Upregulation of NEK2 is associated with drug resistance in ovarian cancer, Oncol. Rep., 31 (2014), 745-754. doi: 10.3892/or.2013.2910
    [49] C. L. Au-Yeung, T. L. Yeung, A. Achreja, H. Zhao, K. P. Yip, S. Y. Kwan, et al., ITLN1 modulates invasive potential and metabolic reprogramming of ovarian cancer cells in omental microenvironment, Nat. Commun., 11 (2020), 3546. doi: 10.1038/s41467-020-17383-2
    [50] R. Januchowski, P. Zawierucha, M. Andrzejewska, M. Ruciński, M. Zabel, Microarray-based detection and expression analysis of ABC and SLC transporters in drug-resistant ovarian cancer cell lines, Biomed. Pharmacother., 67 (2013), 240-245. doi: 10.1016/j.biopha.2012.11.011
    [51] Y. Wang, F. Shao, L. Chen, ALDH1A2 suppresses epithelial ovarian cancer cell proliferation and migration by downregulating STAT3, OncoTargets Ther., 11 (2018), 599-608. doi: 10.2147/OTT.S145864
    [52] S. R. Rosario, M. D. Long, H. C. Affronti, A. M. Rowsam, K. H. Eng, D. J. Smiraglia, Pan-cancer analysis of transcriptional metabolic dysregulation using The Cancer Genome Atlas, Nat. Commun., 9 (2018), 5330. doi: 10.1038/s41467-018-07232-8
    [53] D. Reimer, S. Sadr, A. Wiedemair, S. Stadlmann, N. Concin, G. Hofstetter, et al. Clinical relevance of E2F family members in ovarian cancer-an evaluation in a training set of 77 patients, Clin. Cancer Res., 13 (2007), 144-151. doi: 10.1158/1078-0432.CCR-06-0780
    [54] L. Zhan, Y. Zhang, W. Wang, E. Song, Y. Fan, B. Wei, E2F1: a promising regulator in ovarian carcinoma, Tumor Biol., 37 (2016), 2823-2831. doi: 10.1007/s13277-015-4770-7
    [55] S. Hein, S. Mahner, C. Kanowski, T. Lö ning, F. Jä nicke, K. Milde-Langosch, Expression of Jun and Fos proteins in ovarian tumors of different malignant potential and in ovarian cancer cell lines, Oncol. Rep., 22 (2009), 177-183.
    [56] W. Yang, H. Cho, H. Y. Shin, J. Y. Chung, E. S. Kang, E. J. Lee, et al., Accumulation of cytoplasmic Cdk1 is associated with cancer growth and survival rate in epithelial ovarian cancer, Oncotarget, 7 (2016), 49481-49497. doi: 10.18632/oncotarget.10373
    [57] S. Wei, J. Liu, Y. Shi, X. Zhang, Y. Yang, Q. Song, Exploration of the sequential gene changes in epithelial ovarian cancer induced by carboplatin via microarray analysis, Mol. Med. Rep., 16 (2017), 3155-3160. doi: 10.3892/mmr.2017.7008
    [58] W. Li, Z. Liu, B. Liang, S. Chen, X. Zhang, X. Tong, et al., Identification of core genes in ovarian cancer by an integrative meta-analysis, J. Ovarian Res., 11 (2018), 94. doi: 10.1186/s13048-018-0467-z
    [59] L. Wu, Z. H. Ling, H. Wang, X. Y. Wang, J. Gui, Upregulation of SCNN1A promotes cell proliferation, migration, and predicts poor prognosis in ovarian cancer through regulating epithelial-mesenchymal transformation, Cancer Biother. Radiopharm., 34 (2019), 642-649.
    [60] C. Chen, S. Chen, M. Luo, H. Yan, L. Pang, C. Zhu, et al., The role of the CDCA gene family in ovarian cancer, Ann. Transl. Med., 8 (2020), 190. doi: 10.21037/atm.2020.01.99
    [61] Y. Li, M. Xiao, F. Guo, The role of Sox6 and Netrin-1 in ovarian cancer cell growth, invasiveness, and angiogenesis, Tumour Biol., 39 (2017), 1010428317705508.
    [62] R. Aguirre-Gamboa, H. Gomez-Rueda, E. Martínez-Ledesma, A. Martínez-Torteya, R. Chacolla-Huaringa, A. Rodriguez-Barrientos, et al., SurvExpress: an online biomarker validation tool and database for cancer gene expression data using survival analysis, PloS One, 8 (2013), e74250.
    [63] R. W. Tothill, A. V. Tinker, J. George, R. Brown, S. B. Fox, S. Lade, et al., Novel molecular subtypes of serous and endometrioid ovarian cancer linked to clinical outcome, Clin. Cancer Res., 14 (2008), 5198-5208. doi: 10.1158/1078-0432.CCR-08-0196
    [64] Y. Liu, Survival correlation of immune response in human cancers, Oncotarget, 10 (2019), 6885-6897. doi: 10.18632/oncotarget.27360
    [65] J. Wang, F. H. C. Cheng, J. Tedrow, W. Chang, C. Zhang, A. K. Mitra, Modulation of immune infiltration of ovarian cancer tumor microenvironment by specific subpopulations of fibroblasts, Cancers, 12 (2020), 3184. doi: 10.3390/cancers12113184
    [66] W. Wang, W. Zou, J. R. Liu, Tumor-infiltrating T cells in epithelial ovarian cancer: predictors of prognosis and biological basis of immunotherapy, Gynecol. Oncol., 151 (2018), 1-3. doi: 10.1016/j.ygyno.2018.09.005
    [67] S. Lieber, S. Reinartz, H. Raifer, F. Finkernagel, T. Dreyer, H. Bronger, et al., Prognosis of ovarian cancer is associated with effector memory CD8+ T cell accumulation in ascites, CXCL9 levels and activation-triggered signal transduction in T cells, Oncoimmunology, 7 (2018), e1424672.
    [68] S. U. Kumar, D. T. Kumar, R. Siva, C. G. P. Doss, H. Zayed, Integrative bioinformatics approaches to map potential novel genes and pathways involved in ovarian cancer, Front. Bioeng. Biotechnol., 7 (2019), 391. doi: 10.3389/fbioe.2019.00391
    [69] S. Mishra, M. I. Shah, S. Udhaya Kumar, D. Thirumal Kumar, C. Gopalakrishnan, A. M. Al-Subaie, et al., Chapter Eleven - Network analysis of transcriptomics data for the prediction and prioritization of membrane-associated biomarkers for idiopathic pulmonary fibrosis (IPF) by bioinformatics approach, 123 (2020), 241-273.
    [70] H. Yan, G. Zheng, J. Qu, Y. Liu, X. Huang, E. Zhang, et al., Identification of key candidate genes and pathways in multiple myeloma by integrated bioinformatics analysis, J. Cell. Physiol., 234 (2019), 23785-23797. doi: 10.1002/jcp.28947
    [71] S. Udhaya Kumar, D. Thirumal Kumar, R. Siva, C. George Priya Doss, S. Younes, N. Younes, et al., Dysregulation of signaling pathways due to differentially expressed genes from the B-cell transcriptomes of systemic lupus erythematosus patients-a bioinformatics approach, Front. Bioeng. Biotechnol., 8 (2020), 276. doi: 10.3389/fbioe.2020.00276
    [72] B. Charbonneau, E. L. Goode, K. R. Kalli, K. L. Knutson, M. S. DeRycke, The immune system in the pathogenesis of ovarian cancer, Crit. Rev. Immunol., 33 (2013), 137-164. doi: 10.1615/CritRevImmunol.2013006813
    [73] R. Nameki, H. Chang, J. Reddy, R. I. Corona, K. Lawrenson, Transcription factors in epithelial ovarian cancer: histotype-specific drivers and novel therapeutic targets, Pharmacol. Ther., 220 (2021), 107722. doi: 10.1016/j.pharmthera.2020.107722
    [74] H. Ohtani, Focus on TILs: prognostic significance of tumor infiltrating lymphocytes in human colorectal cancer, Cancer Immun., 7 (2007), 4.
    [75] D. Yang, Y. He, B. Wu, Y. Deng, N. Wang, M. Li, et al., Integrated bioinformatics analysis for the screening of hub genes and therapeutic drugs in ovarian cancer, J. Ovarian Res., 13 (2020), 10. doi: 10.1186/s13048-020-0613-2
    [76] H. Feng, Z. Y. Gu, Q. Li, Q. H. Liu, X. Y. Yang, J. J. Zhang, Identification of significant genes with poor prognosis in ovarian cancer via bioinformatical analysis, J. Ovarian Res., 12 (2019), 35. doi: 10.1186/s13048-019-0508-2
    [77] X. Yang, S. Zhu, L. Li, L. Zhang, S. Xian, Y. Wang, et al., Identification of differentially expressed genes and signaling pathways in ovarian cancer by integrated bioinformatics analysis, OncoTargets Ther., 11 (2018), 1457-1474. doi: 10.2147/OTT.S152238
  • mbe-18-05-324-supplementary.xlsx
  • 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
  • © 2021 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(5041) PDF downloads(179) Cited by(4)

Figures and Tables

Figures(7)  /  Tables(4)

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog