Increased content of monounsaturated fatty acids with low melting points in cattle milk and meat known to strongly affect their taste and nutritional value. Steroyl-CoA-desaturase (SCD) is a key enzyme in the biosynthesis of monounsaturated fatty acids. The gene encoding SCD in cattle is located on chromosome 26 and consists of 6 exons and 5 introns and is expressed in the cells of cow udder, adipose, muscles and other tissues. To date, several SNPs in cattle SCD gene have been found. For SNP rs41255693 (g.10329T > C), a correlation was shown between the content of monounsaturated fatty acids in cattle milk/meat and milk productivity. This article describes the development of a real-time PCR assay for genotyping SNP rs41255693 in cattle SCD gene with allele-specific TaqMan probes. The use of this assay allows to significantly reduce the time needed for cattle genotyping compared to the commonly used PCR-RFLP method.
1.
Introduction
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.
2.
Preliminaries
This section provides definitions, lemmas and methods that will be used in latter sections.
2.1. Concept of fuzzy set
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)):x∈X}. 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 x∈X 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)=¯{x∈R,μ˜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
Therefore, the α-cut of triangular fuzzy number is a bounded and closed internal [AL(α),AR(α)], where AL(α)=inf{x:μ˜A(x)≥α}=a1+α(a2−a1) and AR(α)=sup{x:μ˜A(x)≥α}=a3−α(a3−a2).
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 ˜A∈E, 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,˜B∈E, 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αL−BαR,AαR−Bα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) [limk→∞AαkL,limk→∞Aα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 ˜A∈E, then the conditions (1) and (2) hold true.
2.2. Weighted sum method
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
where wi>0 and ∑liwi=1 are satisfied.
3.
Model foundation
Motivated by the predator-prey system (3) in Wang et al. [5] as follows
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
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(1−m)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
By the method of weighted sum, the above system translates into
where w1 and w2 express as weight satisfying w1+w2=1, and w1,w2≥0. The simplified form of system (3.4) is given by
where
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].
4.
Positivity and boundedness of system
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
which ensures the solution of system (3.5) is positive.
On the other side, the first differential equation of system (3.5) provides
which implies lim supt→∞x(t)≤A1KA2≡k1. From the second differential equation of system (3.5), one has
which can be written as
It is easy to see lim supt→∞y(t)≤k1B3(1−m)B2. Therefore, the solution of system (3.5) is bounded.
5.
Existence and stability of equilibria
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(A1−q1E1)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
exists if B3(1−m)x∗>B1+q2E2 and A1B2+A3(1−m)(B1+q2E2)>B2q1E1 are satisfied.
5.1. Local stability
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
or saddle point when
(3) Interior equilibrium S∗ is stable node or focus when
Proof. The following Jacobian matrix of system (3.5) is expressed as
(1) The Jacobian matrix M0 at S0 is represented as
According to the characteristic equation about the above matrix
one has λ1=A1−q1E1, λ2=−B1−q2E2, hence A0 is stable node when A1<q1E1 or saddle point when A1>q1E1.
(2) The Jacobian matrix M1 at S1 is depicted as
In accordance with the characteristic equation of the above matrix
we have λ1=−(A1−q1E1) and λ2=−(B1+q2E2)+KB3(1−m)(A1−q1E1)A2. Obviously, S1 exists if A1>q1E1 holds, hence S1 is a stable node when A1>q1E1 and A2(B1+q2E2)>KB3(1−m)(A1−q1E1) or saddle point when A1>q1E1 and A2(B1+q2E2)<KB3(1−m)(A1−q1E1).
(3) The Jacobian matrix M∗ at S∗ is written as
The characteristic equation is expressed as
It is easy to see the sum of the roots of Eq (5.6) is negative, i.e., λ1+λ2=−(B2y∗+A2x∗K)<0 and the product of the roots of Eq (5.6) is positive, i.e., λ1λ2=A3B3(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 B3(1−m)x∗>B1+q2E2 and A1B2+A3(1−m)(B1+q2E2)>B2q1E1 are satisfied.
5.2. Global stability
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
where l>0 is a constant. Differentiating both sides of the above equation in regard to time t, one has
Choose l=A3B3, it yields to
Since dVdt in some neighborhood of (x∗,y∗) is negative semidefinite, the equilibrium (x∗,y∗) is globally asymptotically stable.
6.
Bionomic equilibria
As we have discussed in Section 5, the biological equilibrium yields to dxdt≡0 and dydt≡0. 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
where N1=(p1q1x−c1)E1,N2=(p2q2y−c2)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
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
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
which exists if satisfying A1p2q2>A3(1−m)c2. Substituting Eq (6.4) into the third equation of Eq (6.3) gives
which is positive provided that
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
that (x∞,y∞,E1∞,E2∞)=(c1p1q1,0,A1p1q1K−A2c1p1q21K,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
satisfying B3c1(1−m)>B1p1q1. Substituting Eq (6.8) into the second equation of Eq (6.7) gives
which is positive provided that
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
which exist if satisfying the following conditions
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(1−m)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∞=A1p1q1K−A2c1p1q21K.
(3) The boundary bionomic equilibrium (x∞,y∞,E1∞,0) exists if B3c1(1−m)>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).
7.
Fuzzy optimal harvesting
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
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 0≤Ei(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
Again, by the method of weighted sum, it yields to
where w′1,w′2≥0 are two weights such that w′1+w′2=1. The Hamiltonian can be written as
where λi=λi(t)(i=1,2) are the adjoint variables. Applying Pontryagin's maximum principle [60], we get the adjoint equations
It follows from Eqs (7.4) and (7.5) that
Substituting the interior equilibrium into Eq (7.6) yields
Now, we neglect λ1 in Eq (7.7) to get a reduced differential equation for λ2
where
The solution of Eq (7.8) as follows
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
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
Similarly, we get
where
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
that is,
We substitute λ1 and λ2 in Eq (7.11) and Eq (7.10) into Eq (7.14) and obtain
where
Equation (7.15) is equivalent to the following equations
To eliminate e−(δR−δL)t in Eq (7.17), we consider the following two cases:
Case Ⅰ If PL≠c1w′1 and QL≠c2w′1, then Eq (7.17) is equal to
Case Ⅱ If PL=c1w′1 or QL=c2w′1, we have
Motivated by the method in [25], firstly transforming Eq (7.14) into the following form
Differentiating λi(i=1,2) in Eq (7.19) with respect to t, it derives that
Substituting λi(i=1,2) in Eq (7.19) into Eq (7.6) yields
Combined Eqs (7.20) and (7.21) with E1 and E2 omitted, one has
Consider N1=(p1q1x−c1)E1>0 and N2=(p2q2y−c2)E2>0, we divide both sides of the first equation by that of the second equation of Eq (7.22)
Besides, the values of E1 and E2 at the interior equilibrium turn to
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˜δ.
8.
Numerical simulations
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.
Tables 1–3 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 1–3 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.
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 4–6, respectively. Meanwhile, the equilibria are stable under different initial conditions.
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.
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.
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.
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 w′1,w′2,α′ are appeared in Tables 5–9, respectively. For different combinations w1,w2,α and w′1,w′2,α′, we see the existence of optimal equilibrium.
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 w′1; when α′=0.3, x˜δ decreases and then increases, y˜δ is invariant, E1˜δ decreases and then increases, and E2˜δ increases and then decreases with increasing w′1; 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 w′1; when α′=0.9, x˜δ and y˜δ are invariant, E1˜δ decreases and then increases, and E2˜δ is invariant with increasing w′1. Besides, when fixed w′1 and w′2, 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 w′1; when fixed w′1 and w′2, 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 w′1; when fixed w′1 and w′2, 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 w′1; when fixed w′1 and w′2, 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 w′1; when α′=0.3 or α′=0.6, we can see that x˜δ increases and then decreases, y˜δ, E1˜δ and E2˜δ are invariant with increasing w′1; when α′=0.9, x˜δ, y˜δ and E1˜δ are invariant, E2˜δ decreases and then increases with increasing w′1. When fixed w′1 and w′2, x˜δ decreases, y˜δ is invariant, E1˜δ is almost invariant, and E2˜δ increases with the increase of α′.
From Tables 5–7, it is easy to see that when fixed w1,w2 and w′1,w′2,α′, 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 7–9, we observed that when fixed α and w′1,w′2,α′, 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.
9.
Conclusions
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 4–6 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 5–9. We obtain optimal equilibrium and optimal harvesting effort as respect to w′1,w′2 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.
Acknowledgements
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).
Conflict of interest
The authors have no conflict of interest.