1.
Introduction
In population ecology, understanding how predators and primary producers influence nutrient flow relative to each other is important. Ecosystem interactions and predator-prey relationships are governed by predation and the delivery of resource processes. The identification of ecological factors that can alter or control dynamic behavior requires theoretical and experimental research. One way to study these questions is by means of experimental control, and another useful way is via mathematical modeling as well as computer simulations. Over decades of theoretical ecology and biomathematics development, mathematical modeling has become an indispensable tool for scientists in related fields to study ecosystems. Since Lotka [1] and Volterra [2], as cornerstones of theoretical ecology, published the first study of predator-prey dynamics, any species in nature can be a predator or prey, and due to its prevalence, it has become one of the most popular topics for researchers to study [3,4,5]. Besides, because biological resources are renewable and have the most unique development mechanisms, the over-utilization of biological resources and the destruction of the environment by humans will directly affect the balance of the ecosystem. Maintaining ecological balance and meeting humans material needs have attracted the most attention from researchers focused on the scientific management of renewable resource development [6,7,8].
Shelter serves as a defense strategy. It refers broadly to a series of behaviors by prey to avoid predators in order to increase their survival rate. The concept of sanctuary was first developed by Maynard-Smith [9] and Gause et al. [10], and its popularity has been very high, garnering widespread attention from many scholars [11,12,13,14,15]. Sih et al.[16] investigated the effects of prey refuge in a three-species model and concluded that the system's stability is strongly related to the refuge. Also, similar findings can be displayed in [17,18,19,20,21,22]. The two modes of refuge analyzed by Gonzalez-Olivares et al. [17] have diverse stability domains in terms of the parameter space. Qi et al.[21] ensure the stability of the system by varying the strength of the refuge.
Through reviewing a large amount of literature, we begin to consider [23,24] as a basis for the two prey and one predator species that will be modeled in this article. We assume that at a certain time t, the populations of the two prey and one predator are x1(t), x2(t), and y(t), respectively. Based on the above, we construct the following model:
The significance of the full parameters is annotated in Table 1.
Most species in nature, including humans, are influenced by fear. Fear may cause an abnormal state and behavior to arise. As usual, prey have an innate fear of predators. The ecology of fear is related to combining the optimal behavior of prey and predators with their population densities [25,26]. In view of reality, it is a fact that prey fear predators, which is seen as a psychological effect that can have a lasting impact on prey populations. This psychological influence is often easy to overlook, but it is necessary to consider it in the context of practical ecology [27]. Wang et al. [28] first considered the effect of the fear factor on the model and first proposed the fear of prey F(k,y). Afterwards, some researchers have investigated the effects of the fear effect and predator interferences in some three-dimensional systems as well as explored the generation of Hopf bifurcation conditions in the presence of a fear parameter as a bifurcation parameter [29,30,31,32]. Zanette et al. [33] observed that prey will reduce reproducing because of fear of being killed by predators, thus decreasing the risk of being killed after giving birth, which also leads directly to a decline in prey birth rates. According to the above discussion, our paper considers the different fears ki caused by predators for the two prey species.
In reality, when prey feel the crisis of being hunted, they will reproduce less and increase their survival rate. These conditions about the fear factor F(ki,y) (i=1,2) are listed as follows:
1) F(0,y)=1: prey production does not decrease when the prey does not fear the predator;
2) F(ki,0)=1: even though the prey will develop a fear of predators and there will be no predators, prey production will still not decline;
3) limki→∞F(ki,y)=0: when the prey's fear of the predator is very high, this will result in the prey production tending to zero;
4) limy→∞F(ki,y)=0: prey have a fear of predators, and when predator numbers are too large, this can also lead to prey production tending to zero;
5) ∂F(ki,y)∂ki<0: the greater the prey's fear of predators, the less productive it will be;
6) ∂F(ki,y)∂y<0: predators are inversely proportional to their prey.
For ease of analysis, we draw on Wang et al. [28] to consider the fear effect:
obviously, F(ki,y) (i=1,2) in (1.2) satisfies conditions 1)–6). Based on the above conditions, this study will consider the effect of fear on system (1.1) to obtain system (1.3).
Notably, most biological parameters in much of the literature are fixed constants. However, in reality, the survival of species is full of unknowns, and all data are not always constant, which can lead to deviations from the ideal model with fixed parameters. In order to make the model more relevant and the results more accurate, we cannot just consider fixed parameters. Therefore, to make the study more convincing, it is necessary to target imprecise parameters. Professor Zadeh [34], who first proposed the fuzzy set theory, also argued that the application of fuzzy differential equations is a more accurate method for modeling biological dynamics in the absence of accurate data conditions [35]. Moreover, the first introduction of the idea of fuzzy derivatives came from Chang and Zadeh [36]. Further, Kaleva [37] studied the generalized fuzzy derivatives based on Hukuhara differentiability, the Zadeh extension principle, and the strong generalized differentiability concept. Bede et al. [38] employed the notion of strongly generalized differentiability to investigate fuzzy differential equations. Khastan and Nieto [39] solved the margin problem for fuzzy differential equations in their article. Motivated by the method of Pal [13] and Wang [23], we assume that the imprecise parameters ~r1, ~r2, ~a1, ~a2, ~c1, ~c2, ~e1, ~e1 and ˜d represent all triangular fuzzy numbers (the relevant theories of fuzzy sets are detailed in Appendix A), then the system (1.3) can be written as
we cut these imprecise parameters ~r1, ~r2, ~a1, ~a2, ~c1, ~c2, ~e1, ~e1, and ˜d by using α-level. System (1.4) can be expressed as follows:
Introducing weighted sum, we change (1.5) to (1.6)
where w1 and w2 are satisfied with w1+w2=1, and w1,w2≥0. Simplifying the system (1.6), we obtain
where
The rest of the paper is shown below: In Section 2, we first prove the nonnegativity and boundedness of the system (1.7). Sections 3 and 4 discuss all possible equilibria and give conditions for the local asymptotic stability and global asymptotic stability of the equilibria. Immediately after that, in Section 5, we analyze the Hopf bifurcation by using the normal form theory. In Section 6, we numerically simulate the theoretical results of Sections 4 and 5. Finally, the article ends with detailed conclusions.
2.
Nonnegativity and boundedness
In this section, we give the following theorem to ensure the boundedness and nonnegativity of the solutions of the system (1.7).
Theorem 2.1. Provided that the initial values x1(0)>0, x2(0)>0, and y(0)>0, all solutions of system (1.7) are nonnegative.
Proof. It is not difficult to find that the right half of the system (1.7) fulfills the local Lipschitzian condition. Integrating both sides of the system (1.7) at the same time yields
If the solution curve starts at any internal point of R3+={(x1(t),x2(t),y(t))∈R3:x1(t)≥0, x2(t)≥0, y(t)≥0}, then x1(t), x2(t), and y(t) will always be nonnegative. □
Theorem 2.2. Assume that the initial values x1(0), x2(0), and y(0) are all greater than zero. The feasible region Ω is a positive invariant set of the system (1.7) defined by
where μ=min{q1E1,q2E2,C3+q3E3}.
Proof. Define a function
After taking the derivative on both sides of (2.2), we obtain
Furthermore, we can obtain
where μ=min{q1E1,q2E2,C3+q3E3}. Let ϕ1=A21K14A2,ϕ2=B21K24B2, ϕ≡C1A4ϕ1+C2B4ϕ2, we have
Therefore, it can be deduced that
where N is a positive constant. Then we can further obtain
which indicates that the feasible domain Ω is a positive invariant set. □
3.
Existence of biological equilibria
In this section, we discuss the existence of all equilibria in the system (1.7). All equilibria for system (1.7) are provided by
(1) Trivial equilibrium P1=(0,0,0).
(2) Axial equilibrium P2=(xς1,0,0) exists if A1>q1E1, where xς1=K1(A1−q1E1)A2.
(3) Axial equilibrium P3=(0,xΥ2,0) exists if B1>q2E2, where xΥ2=K2(B1−q2E2)B2.
(4) Axial equilibrium P4=(xΨ1,0,yΨ) exists if A1≥A2xΨ1K1+q1E1 and √Δ1>A4(1−m1)+k1q1E1, where
(5) Axial equilibrium P5=(0,xν2,yν) exists if B1≥B2xν2K2+q2E2 and √Δ2>B4(1−m2)+k2q2E2, where
(6) Axial equilibrium P6=(xι1,xι2,0) exists if B1>B3xι1+q2E2, K1K2A3B3>A2B2 and A3B1K2+B2q1E1>A3K2q2E2+A1B2, where
(7) Internal equilibrium P7=(x∗1,x∗2,y∗) exists, and its value will be given in the proof of Theorem 3.1.
Theorem 3.1. When g3>0 and g4g5<0 are met, there is an internal equilibrium P7.
Proof. We derive that from the second equation of the system (1.7)
where
It follows from the Descartes law of signs that Eq (3.4) has one and only one solution y∗ greater than zero if and only if g3>0, i.e., B1>x2K2B2+B3x1+q2E2. Substituting y∗ into the algebra expression on the right side of the first equation of the system (1.7) equals zero; furthermore, we obtain
Introduce (3.5) into the right side of the third equation of the system (1.7), which meets zero, it simplifies to obtain
where
Reusing the Descartes law of signs, we can assert that there exists at least one positive solution x∗1 of Eq (3.6) if and only if g4g5<0. And then we can deduce that
then the interior equilibrium P7(x∗1,x∗2,y∗) exists.□
4.
Stability analysis
In this section, the Jocabian matrix will be used to prove the local stability of all equilibria. Moreover, we prove the global stability of the internal equilibrium P7 by constructing a Lyapunov function.
4.1. Local stability
The Jocabian matrix for system (1.7) is given below:
where
Through simple calculation, we directly draw the conclusion that trivial and axial equilibria are locally asymptotically stable:
(1) P1(0,0,0) is locally asymptotically stable if
(2) P2(xς1,0,0) is locally asymptotically stable if
(3) P3(0,xΥ2,0) is locally asymptotically stable if
(4) P4(xΨ1,0,yΨ) is locally asymptotically stable if
(5) P5(0,xν2,yν) is locally asymptotically stable if
(6) P6(xι1,xι2,0) is locally asymptotically stable if
We draw the conclusion that the internal equilibrium P7(x∗1,x∗2,y∗) is locally asymptotically stable from the proof of Theorem 4.1.
Theorem 4.1. The internal equilibrium P7 is locally asymptotically stable if it exists and the following conditions are fulfilled:
Proof. The Jocabian matrix of system (1.7) at (x∗1,x∗2,y∗) is
where
Therefore, the characteristic equation at P7 can be expressed as
where
The Routh-Hurwitz criterion shows that the internal equilibrium P7 is locally asymptotically stable; the following conditions need to be met: ψ1>0, ψ1ψ2>0, and ψ3>0.□
4.2. Global stability
This subsection studies the global asymptotic stability of interior equilibrium P7.
Theorem 4.2. If condition 4Γ1Γ2l1l2A2B2(1+k1y)(1+k2y) > (l1A3+l2B3)2 (i.e. 4Γ1Γ2l1l2(w1rα1R+w2rα1L)(w1rα2R+w2rα2L)(1+k1y)(1+k2y) > (l1(w1aα1R+w2aα1L)+l2(w1aα2R+w2aα2L))2) holds, then P7 is globally asymptotically stable.
Proof. We construct a Lyapunov function:
Obviously, xi−x∗i−x∗iln(xix∗i)≥0 (i=1,2) and y−y∗−y∗ln(yy∗)≥0, thus V≥0. Taking the derivative of V(x1,x2,y) over t, one has
where
To simplify the calculation, let
We obtain
where
Therefore, dVdt<0 if and only if 4Γ1Γ2l1l2A2B2(1+k1y)(1+k2y)>(l1A3+l2B3)2. □
5.
Hopf bifurcation
In this section, we will use the normal form theory introduced by Hassard et al.[40] and the central manifold theory [41] to study the Hopf bifurcation of the system (1.7). When the system (1.7) undergoes Hopf bifurcation, the corresponding characteristic equation must have a pair of conjugate pure imaginary roots, that is,
Consider the parameter k1 as a bifurcation parameter. When the value of parameter k1 changes near the critical point kΞ1 of Hopf bifurcation, the pure imaginary roots ±iω will become a complex eigenvalue η=ρ+i˜ω. Substituting η=ρ+i˜ω into Eq (4.12), we need to separate the imaginary and real parts to get
By simplifying Eqs (5.2) and (5.3), we obtain
at k1=kΞ1, taking the derivative of Eq (5.4) over k1 yields
If it satisfies dρdk1|k1=kΞ1≠0, the system (1.7) will generate Hopf bifurcation, which indicates that when parameter k1 crosses the bifurcation critical point kΞ1, the population state evolves from stable equilibrium to periodic oscillation over time.
When the system (1.7) undergoes Hopf bifurcation at k1=kΞ1, the final decision condition is also met. Considering that the characteristic roots of Eq (4.12) are η1,2=±iω and η3=−ψ1, in order to obtain this condition, we introduce
Substituting (5.6) into the system (1.7) and separating the linear and nonlinear parts, it can be obtained that
where
where O((|z1|+|z2|+|z3|)4) is a fourth-order polynomial function about variables (|z1|,|z2|,|z3|), while tj1j2j3, nj1j2j3, and lj1j2j3 can be obtained through calculation:
By introducing a reversible transformation
which q21q32−q22q31+q22q33−q23q32≠0, the expression for the coefficient qij(i=2,3;j=1,2,3) is
Li(i=1,2,⋯,9) is defined in (4.11). So the standard type of system (5.7) can be written as
where
In Eq (5.13), the coefficients of polynomial ˜F(y1,y2,y3) are ˜tj1j2j3, ˜nj1j2j3 and ˜lj1j2j3.
Based on the central manifold theory, use the center manifold Wc(0,0,0) existing at the origin to reduce the dimension of the system (5.12), that is
where
Through calculation, we can obtain
Correspondingly, the dynamic properties of the system are limited to the central flow Wc(0,0,0), and in conjunction with Eq (5.14), system (5.12) can be simplified as
where
in the formula, we have
We introduce the partial derivative sign
where subscripts y1 and y2 indicate partial derivatives for the first and second variable, respectively. Based on Eq (5.18), it can be obtained that Uy1=0, Uy2≠0, Ny1≠0, Ny2=0, and Uy2Ny1≠0. In addition, it ensures that the system (5.18) has pure virtual feature roots ±i√|Uy2Ny1|. Thus, it can be determined that system (1.7) produces Hopf bifurcation; the direction of the bifurcation is determined by the following equation:
Theorem 5.1. If dρdk1|k1=kΞ1≠0, then system (1.7) will generate Hopf bifurcation at interior equilibrium P7. In addition, when dρdk1|k1=kΞ1<0, if QkΞ1<0 and 0<k1−kΞ1≪1, then system (1.7) will generate supercritical Hopf bifurcation and form a stable periodic orbit, or if QkΞ1>0 and 0<k1−kΞ1≪1, then system (1.7) will generate subcritical Hopf bifurcation and form a stable periodic orbit.
6.
Numerical simulations
In this section, we first discussed equilibria P1 to P7 of system (1.7) with distinct values of α, w1, and w2. Consider the parameter values as follows: ~r1=(2.8,3,3.2), ~r2=(2.8,3,3.2), ~c1=(0.1,0.2,0.3), ~c2=(0.5,0.6,0.7), ~a1=(0.1,0.2,0.3), ~a2=(0.2,0.3,0.4), ~e1=(0.2,0.3,0.4), ~e2=(0.3,0.4,0.5), and ˜d=(0.1,0.2,0.3). Tables 2–8 showed that the trivial equilibrium P1 retained constant at (0, 0, 0), the values of prey x1, prey x2, and predator y always maintained at 0; the values of prey x1 in P2 and prey x2 in P3 severally decreased with increasing w1 under the same α; the values of prey x1 and predator y in P4 increased with increasing w1, and for P5 the value of prey x2 and predator y rose with growing w1; the values of prey x1 and x2 in P6 decreased with growing w1; and for the same α, considering interior equilibrium P7, the values of prey x1, prey x2, and predator y decreased with growing w1.
Considering four sets of different initial values, it could be seen from Figure 1 that different orbits eventually converged to the same value, which concluded that the interior equilibrium of the system (1.7) fulfills the character of globally asymptotical stability. Figure 2 plotted the bifurcation graph of system (1.7) with the horizontal coordinates k1, and the Hopf bifurcation of the system occurred with k1 taking values in the range of 0.01≤k1≤0.7. When 0.01≤k1<0.384, the system oscillates periodically, while it maintains a stable steady-state when 0.384<k1≤0.7. Therefore, based on Figure 2, it could be concluded that the fear of prey x1 for predator y affected the stability of the system. We further observed that as k1 increased, the prey x1 density continued to decrease while the predator y density kept increasing. Thus, the result also suggested that greater fear of predators had a negative impact on prey populations while having a positive impact on predator populations. Correspondingly, Figures 3 and 4 showed the waveform plots and phase diagram at k1=0.1 and k1=0.7, respectively.
In addition, Figure 5 also plots the bifurcation graph with changing m1. As can be seen in Figure 5, m1 took values from 0.3 to 1, in which the system also underwent a Hopf bifurcation. When the value m1 ranged from 0.3 to 0.657, the system (1.7) was stable; nevertheless, it would become unstable at 0.657<m1≤1. Correspondingly, Figures 6 and 7 showed the waveform plots and phase diagram at m1=0.6 and m1=0.9, respectively.
Further, we find an interesting dynamic phenomenon through some numerical simulations. System (1.7) appears as a chaotic phenomenon, as shown in Figure 8.
7.
Conclusions
In this work, we develop a model of one-predator and two-prey interactions in a fuzzy environment, examine the effects of fear and prey refuge on the system, and provide insight into the dynamic complexity. The proofs of the theoretical parts of this paper are based on system (1.7). It has been proven that all equilibria in system (1.7) are locally asymptotically stable, and interior equilibrium P7 is also globally asymptotically stable. We have been further concerned about the appearance and direction of Hopf bifurcation. With the support of theoretical research, our numerical simulations have been able to display a wealth of charts and graphs.
First of all, different equilibria are displayed from Tables 2–8 with different α,w1,w2, respectively. Throughout Figure 1, we have verified the global asymptotical stability of interior equilibrium P7, and find that the system is from unstable to stable with the increase of fear k1, which demonstrates that the fear effect may be an important factor influencing the stability of the system (see Figures 2–4). Furthermore, it has also been observed that an increase in prey refuge m1 leads to oscillatory phenomena (see Figures 5–7). Finally, through studying the Hopf bifurcation, we have discovered some interesting biological phenomena, namely that system (1.7) appears to be in a chaotic state (see Figure 8).
Author contributions
Xuyang Cao: Conceptualization, Investigation, Methodology, Validation, Writing-original draft, Formal analysis, Software; Qinglong Wang: Conceptualization, Methodology, Formal analysis, Writing-review and editing, Supervision; Jie Liu: Validation, Visualization, Data curation.
Use of AI tools declaration
The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.
Acknowledgments
The authors thank the editor and referees for their careful reading and valuable comments.
The work is supported by the Natural Science Foundation of Hubei Province (No. 2023AFB1095) and 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 Teaching Research Project of Education Department of Hubei Province (No. 2022367) and the Graduate Education Innovation Project of Hubei Minzu University (Nos. MYK2024071, MYK2023042).
Conflict of interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Appendix A
Definition 1. [34] Fuzzy set: A fuzzy set ˜ℏ in a universe of discourse S is denoted by the set of pairs
where the mapping μ˜ℏ:S→[0,1] is the membership function of the fuzzy set ˜ℏ and μ˜ℏ is the membership value or degree of membership of s∈S in the fuzzy set ˜ℏ.
Definition 2. [42] α-cut of fuzzy set: For any α∈(0,1], the α-cut of fuzzy set ˜ℏ defined by ℏα={s:μ˜ℏ(s))≥α} is a crisp set. For α=0 the support of ˜ℏ is defined as ℏ0=Supp(˜ℏ)=¯{s∈R,μ˜ℏ(s)>0}.
Definition 3. [43] Fuzzy number: A fuzzy number satisfying the property S=R is called a convex fuzzy set.
Definition 4. [44] Triangular fuzzy number: A triangular fuzzy number (TFN) ˜ℏ≡(b1,b2,b3) represent fuzzy set of the real line R satisfying the property that the membership function μ˜ℏ:R→[0,1] can be espressed by
Hence, the α-cut of triangular fuzzy number meets boundedness and encapsulation on [ℏL(α),ℏR(α)], in which ℏL(α)=infs:μ˜ℏ(s)≥α=b1+α(b2−b1) and ℏR(α)=sup{s:μ˜ℏ(s)≥α}=b3+α(b3−b2).
Lemma 1. [45] In weighted sum method, wj stands for the weight of jth objective. wjgj represent a utility function for jth objective, and the total utility function π is represented by
where wj>0 and ∑ljwj=1 are satisfied.