Skip to main content

ORIGINAL RESEARCH article

Front. Appl. Math. Stat., 06 January 2023
Sec. Mathematical Biology
This article is part of the Research Topic Justified Modeling Frameworks and Novel Interpretations of Ecological and Epidemiological Systems View all 11 articles

Bifurcation analysis of a predator–prey model involving age structure, intraspecific competition, Michaelis–Menten type harvesting, and memory effect

  • Biomathematics Research Group, Department of Mathematics, Universitas Negeri Gorontalo, Gorontalo, Indonesia

The complexity of the dynamical behaviors of interaction between prey and its predator is studied. The prey and predator relationship involves the age structure and intraspecific competition on predators and the nonlinear harvesting of prey following the Michaelis–Menten type term. Some biological validities are shown for the constructed model such as the existence and uniqueness as well as the non-negativity and boundedness of solutions. Three equilibrium points, namely the origin, axial, and interior points, are found including their global dynamics by employing the Lyapunov function along with the generalized Lassale invariant principle. The changes in dynamical behaviors driven by the harvesting and the memory effect are exhibited, including transcritical, saddle-node, backward, and Hopf bifurcations. The appearance of these interesting phenomena is strengthened by giving numerical simulations consisting of bifurcation diagrams, phase portraits, and their time series.

1. Introduction

Since Lotka and Volterra introduced the classical predator–prey model, theoretical studies of predation without age structure have attracted the attention of many authors, for example Deng et al. [1], Huang et al. [2], Tahara et al. [3], and Zeng et al. [4]. However, in nature, many species of plants and animals could have life histories that can simply be partitioned into two age stages: immature and mature stages. In each stage, individuals of species have identical biological characteristics, such as the ability to reproduce, motile, ingest food, and survive [5]. In particular, there are amphibians, insects, birds, and mammals with life cycles that can last from only several days or weeks to more than a century. For this reason, some researchers have developed the predator–prey model by incorporating age structure either in prey or/and predator population with other factors that also influence the dynamics of the predator–prey model, mainly restricted to the classical integer-order, stochastic, or delay equations [613].

In Wang and Chen [14] considered the predator–prey model with age structure for the predator population using time delays. If we ignore the effect of time delay, the model can be written as follows:

dxdt= rx(1-xK)-mxz,dydt= nxz-βy-δ1y,dzdt= βy-δ2z.    (1)

Here x(t), y(t), and z(t) represent the population densities of prey, immature predator, and mature predator, at time t, respectively. Model (Equation 1) assumes that the prey grows logistically with r as the intrinsic growth rate, K is the carrying capacity; m is the linear Holling type I functional response, n is the conversion rate with which captured prey are converted to new immature predator, β is the maturity rate of the predator, δ1 and δ2 are the death rate of the immature and mature predators, respectively. It is also assumed that only the mature predator can feed the prey through the term mxz. If we do not consider the age structure of the predator population, then model (Equation 1) is reduced to the classical Lotka–Volterra model for which the positive equilibrium or the boundary equilibrium of this model is globally asymptotically stable. This means that the model has no periodic solution. On the other hand, Wang and Chen [14] prove that in the model (Equation 1), there exists an orbitally asymptotically stable periodic solution around the interior equilibrium point which suggests that the age structure can cause periodic oscillation of populations.

From the point of view of human needs, harvesting of populations generally occurs in wildlife, forestry, and fisheries management. When harvesting is integrated into the predator–prey model, there are three types of harvesting, namely constant harvesting [15], linear harvesting [16], and non-linear harvesting [17]. In this article, we assume that the predator is not a commercial species and there is intraspecific competition among immature predators. Therefore, the predator–prey model with age structure and intraspecific competition in predator (Equation 1), where the prey population is subject to Michaelis–Menten type harvesting, is given by

dxdt= rx(1-xK)-mxz-hxc+x,dydt= nxz-βy-δ1y-ωy2,dzdt= βy-δ2z.    (2)

An example of prey-predator interactions whose biological phenomena are described in the model (Equation 2) can be found in the African wild dog with its prey impala. The African wild dogs are a social structure that lives in packs. For 3–4 weeks, young African wild dogs were in the den with their mother. All adult members of African wild dogs are care for their young ones and provide food for them. The hunting members of the pack will return to the den where they regurgitate meat for the nursing female and young. In some cases, young ones fail to survive because the hunting member does not bring back sufficient food for the young, which leads to intraspecific competition in immature predator [18]. On the other hand, the prey, impala, even though there are no major threats to their survival, poaching has become significantly contributed to the decline in its number [19].

Note that the growth rates of the prey, immature, and mature predator populations in the model (Equation 2) depend only on the local state as the left-hand side is the integer-order derivative. On the other hand, most biological systems have properties where the current state is affected by all of the past states or it is called the memory effect. Therefore, modeling with memory effects can be done by analyzing the system using fractional-order calculus [20, 21]. The operators of the fractional-order derivative have non-local properties to make them more suitable for dynamical systems that have memory influences on their state variables.

After Riemann and Liouville generalized the concept of integer-order calculus to the fractional-order calculus over two decades ago, the studies about the predator–prey models with fractional-order differential equation have gained much attention, for example, Rahmi et al. [22], Owolabi [23], Barman et al. [24], Ghanbari and Djilali [25], Yousef et al. [26], Ghosh et al. [27], and Panigoro et al. [28]. The fractional-order derivatives are defined as an integration that provides the ability to store the whole memory over time, and hence, it could give an exact description of different ecological phenomena. For this reason, the new structure for the model (Equation 2) is given in the following form:

CDtαx(t)= rx(1-xK)-mxz-hxc+x=F1(x,y,z),CDtαy(t)= nxz-βy-δ1y-ωy2=F2(x,y,z),CDtαz(t)= βy-δ2z=F3(x,y,z).    (3)

The existence and local stability of all equilibrium points of the model (Equation 3) are discussed in Panigoro et al. [29]. However, to the best of our knowledge, the global dynamics and bifurcation analysis of the model (Equation 3), to this day, have not been investigated. Here, CDtαf(t) is the standard Caputo derivative for a continuous function f(x) ∈ ℂ([0, +∞), ℝ), which is defined as follows:

 CDtαf(t)=1Γ(1-α)0tf(τ)(t-τ)αdτ,    (4)

where Γ(∗) is the gamma function, t ≥ 0, and 0 < α ≤ 1 is known as the order of the fractional derivative.

Based on the above discussion, we have organized our work in several sections: In Section 3, we develop the existence and uniqueness solution of the model (Equation 3). To check the biologically well-posedness of the model, we establish the non-negativity and boundedness of solutions of the model in Section 3. In Section 4, we derive some sufficient conditions to ensure the global asymptotic stability of each locally asymptotically stable equilibrium point by applying the Lyapunov functions. We then analyze the existing conditions of transcritical, saddle-node, backward, and Hopf bifurcations in Section 5. Some numerical simulations of our obtained results are carried out in Section 6. Finally, the conclusions are given in Section 7.

2. Existence and uniqueness

In this section, we will show that the model (Equation 3) has a unique solution. A similar manner given by Mahata et al. [30] is adopted. We first take the Riemann–Liouville integral (Definition 1 in Yavuz and Sene [31]) on both sides of Equation (3) to achieve the following Volterra-type integral equations.

x(t)-x(0)= 1Γ(α)0tF1(x(τ),y(τ),z(τ))(t-)α-1 dτ,y(t)-y(0)= 1Γ(α)0tF2(x(τ),y(τ),z(τ))(t-τ)α-1 dτ,z(t)-z(0)= 1Γ(α)0tF3(x(τ),y(τ),z(τ))(t-τ)α-1 dτ,    (5)

Now, we will show that the kernels Fi(x, y, z), i = 1, 2, 3, satisfy the Lipschitz condition. For ∥.∥ is the Euclidean norm, we suppose that ∥x(t)∥ ≤ a1, x̄(t)a2, ∥y(t)∥ ≤ a3, ∥ȳ(t)∥ ≤ a4, ∥z(t)∥ ≤ a5, and z̄(t)a6 are bounded functions. For x, x̄, y, ȳ, z, and z̄, we have

      F1(x,y,z)-F1(x̄,y,z)= ||[rx(1-xK)-mxz-hxc+x]-[rx̄(1-x̄K)-mx̄z-hx̄c+x̄]||= ||r(x-x̄)-rK(x+x̄)(x-x̄)-mz(x-x̄)      -ch(x-x̄(c+x)(c+x̄))|| rx-x̄+(a1+a2)rKx-x̄+a5mx-x̄+hcx-x̄= g1x-x̄,      F2(x,y,z)-F2(x,ȳ,z)= [nxz-βy-δ1y-ωy2]-[nxz-βȳ-δ1ȳ-ωȳ2]= -β(y-ȳ)-δ1(y-ȳ)-ω(y+ȳ)(y-ȳ) βy-ȳ+δ1y-ȳ+ω(a3+a4)y-ȳ= g2y-ȳ,      F3(x,y,z)-F3(x,y,z̄)= [βy-δ2z]-[βy-δ2z̄] g3z-z̄,

where g1=r+(a1+a2)rK+a5m+hc, g2 = β + δ1 + ω(a3 + a4), and g3 = δ2. Therefore, we conclude that Fi(x, y, z), i = 1, 2, 3, satisfy the Lipschitz condition. Furthermore, it is clear that if 0 ≤ gi < 1, then Fi(x, y, z) are contractions for i = 1, 2, 3. Therefore, the following theorem is obtained.

Theorem 1. The kernel Fi(x, y, z), i = 1, 2, 3 satisfy the Lipschitz condition and contractions if 0 < gi < 1, i = 1, 2, 3.

Next, Equation (5) can be written as follows:

x(t)= x(0)+1Γ(α)0tF1(x(τ),y(τ),z(τ))(t-τ)α-1 dτ,y(t)= y(0)+1Γ(α)0tF2(x(τ),y(τ),z(τ))(t-τ)α-1 dτ,z(t)= z(0)+1Γ(α)0tF3(x(τ),y(τ),z(τ))(t-τ)α-1 dτ,

Which can be written by the following recursive formula

xn(t)= x(0)+1Γ(α)0tF1(xn-1,y,z)(t-τ)α-1 dτ,yn(t)= y(0)+1Γ(α)0tF2(x,yn-1,z)(t-τ)α-1 dτ,zn(t)= z(0)+1Γ(α)0tF3(x,y,zn-1)(t-τ)α-1 dτ,

with initial conditions x0(t) = x(0), y0(t) = y(0), and z0(t) = z(0). Therefore, we have

φ1n(t)= xn(t)-xn-1(t)=1Γ(α)0t(F1(xn-1,y,z)-F1(xn-2,y,z))(t-τ)α-1 dτ,φ2n(t)= yn(t)-yn-1(t)=1Γ(α)0t(F2(x,yn-1,z)-F2(x,yn-2,z))(t-τ)α-1 dτ,φ3n(t)= zn(t)-zn-1(t)=1Γ(α)0t(F3(x,y,zn-1)-F3(x,y,zn-2))(t-τ)α-1 dτ,    (6)

where xn(t)=j=1nφ1n(t), yn(t)=j=1nφ2n(t), and zn(t)=j=1nφ3n(t). Now, we evaluate the norm of Equation (6). We achieve

φ1n(t)= xn(t)-xn-1(t)                      1Γ(α)0t||(F1(xn-1,y,z)                      -F1(xn-2,y,z))(t-τ)α-1 dτ,φ2n(t)= yn(t)-yn-1(t)                      1Γ(α)0t||(F2(x,yn-1,z)                      -F2(x,yn-2,z))(t-τ)α-1 dτ||,φ3n(t)= zn(t)-zn-1(t)                      1Γ(α)0t||(F3(x,y,zn-1)                      -F3(x,y,zn-2))(t-τ)α-1 dτ||.    (7)

From Theorem 1, we have that the kernel satisfy the Lipschitz condition and hence Equation (7) becomes

xn(t)-xn-1(t)g1Γ(α)0txn-1-xn-2(t-τ)α-1 dτ,yn(t)-yn-1(t)g2Γ(α)0tyn-1-yn-2(t-τ)α-1 dτ,zn(t)-zn-1(t)g3Γ(α)0tzn-1-zn-2(t-τ)α-1 dτ.

The last inequality gives

φ1n(t)g1Γ(α)0tφ1n-1(τ) dτ,φ2n(t)g2Γ(α)0tφ2n-1(τ) dτ,φ3n(t)g3Γ(α)0tφ3n-1(τ) dτ.    (8)

Finally, the existence of a solution is given by the following theorem.

Theorem 2. The solution of model (Equation 3) has a solution under the condition if we have t1 such that (t1giΓ(α+1))<1, i=1,2,3.

Proof. We assume that x(t), y(t), and z(t) are bounded and their kernel Fi, i = 1, 2, 3 satisfy the Lipschitz condition. According to Equation (8), we obtain

φ1n(t) x0(t)[t1g1Γ(α+1)]n,φ2n(t) y0(t)[t1g2Γ(α+1)]n,φ3n(t) z0(t)[t1g3Γ(α+1)]n,    (9)

which represent the existence and continuity of the system. To show that the solution of the model (Equation 3) can be set up from the functions in Equation (9), we assume

x(t)-x(0)=xn(t)-Q1n(t),y(t)-y(0)=yn(t)-Q2n(t),z(t)-z(0)=zn(t)-Q3n(t).    (10)

where Qin, i = 1, 2, 3 are the remaining terms. Furthermore, the given terms would be demonstrated to hold Qi → 0 ∀i = 1, 2, 3. Denoting that

Q1n(t) 1Γ(α)0t(F1(x,y,z)-F1(xn-1,y,z)) dτ                    1Γ(α)0tF1(x,y,z)-F1(xn-1,y,z) dτ                    tg1Γ(α+1)x-xn-1,Q2n(t) 1Γ(α)0t(F2(x,y,z)-F2(x,yn-1,z)) dτ                     1Γ(α)0tF2(x,y,z)-F2(x,yn-1,z) dτ                     tg2Γ(α+1)y-yn-1.Q3n(t) 1Γ(α)0t(F3(x,y,z)-F3(x,y,zn-1)) dτ                     1Γ(α)0tF3(x,y,z)-F3(x,y,zn-1) dτ                     tg3Γ(α+1)z-zn-1.    (11)

By applying a recursive pattern, we acquire

Q1n(t)[tΓ(α+1)]n+1g1nk,Q2n(t)[tΓ(α+1)]n+1g2nk,Q3n(t)[tΓ(α+1)]n+1g3nk.    (12)

At the point t1, we have

Q1n(t)[t1Γ(α+1)]n+1g1nk,Q2n(t)[t1Γ(α+1)]n+1g2nk,Q3n(t)[t1Γ(α+1)]n+1g3nk.    (13)

By considering the results of Theorem 1 and applying n → ∞ to both sides, we have ∥Qi∥ → 0 ∀i = 1, 2, 3.

In the end, we will show that the solution is unique for each initial value by utilizing the contradiction approach. Suppose that there exists another solution of the model (Equation 3), namely x1(t), y1(t), and z1(t). Then, we have

x(t)-x1(t)= 1Γ(α)0t(F1(x,y,z)-F1(xn-1,y,z)) dτ,y(t)-y1(t)= 1Γ(α)0t(F2(x,y,z)-F2(x,yn-1,z)) dτ,z(t)-z1(t)= 1Γ(α)0t(F3(x,y,z)-F3(x,y,zn-1)) dτ.    (14)

Applying the norm on both sides, we achieve

x(t)-x1(t) 1Γ(α)0tF1(x,y,z)-F1(xn-1,y,z) dτ,y(t)-y1(t) 1Γ(α)0tF2(x,y,z)-F2(x,yn-1,z) dτ,z(t)-z1(t) 1Γ(α)0tF3(x,y,z)-F3(x,y,zn-1) dτ.    (15)

By considering Theorem 1, we obtain

x(t)-x1(t) tg1Γ(α+1)x(t)-x1(t),y(t)-y1(t) tg2Γ(α+1)y(t)-y1(t),z(t)-z1(t) tg3Γ(α+1)z(t)-z1(t).

Therefore, the following equations are obtained.

x(t)-x1(t)(1-tg1Γ(α+1)) 0,y(t)-y1(t)(1-tg2Γ(α+1)) 0,z(t)-z1(t)(1-tg3Γ(α+1)) 0.    (16)

As a result, we achieve ∥x(t) − x1(t)∥ = 0, ∥y(t) − y1(t)∥ = 0, and ∥z(t) − z1(t)∥ = 0, which impact x(t) = x1(t), y(t) = y1(t), and z(t) = z1(t). Then, we finally give the following theorem.

Theorem 3. The Caputo fractional-order predator–prey model (Equation 3) has a unique solution.

3. Non-negativity and boundedness

In this section, we will show that for any initial condition is in +3 where

+3:={(x,y,z) : x0, y0, z0, x, y}.

The solution not only exists and is unique but also bounded and always in +3 as t → ∞. Therefore, we have the following two theorems.

Theorem 4. If the initial condition in +3 then both population densities of prey and predator given by model (Equation 3) remain in +3.

Proof. To prove this non-negativity condition, we apply reductio ad absurdum (contradiction method), which is also used in Barman et al. [24] and Maji [32]. We assume that there exists t^>0 such that

{x(t)    >0, when 0t<t^,x(t^)    =0, when t=t^x(t^+)<0, when t=t^+.    (17)

From the first equation in Equation (3) along with Equation (17), we have

 CDtαx(t)|x(t^)=0.    (18)

According to Lemma 3.1 in Barman et al. [33], we get x(t^+)=0 which contradicts with Equation (17) where x(t^+)<0. This means that x(t) ≥ 0 for all t ∈ [0, ∞]. Similarly, we can show that y(t) ≥ 0 and z(t) ≥ 0 for all t ∈ [0, ∞]. In conclusion, we have the non-negative solution for model (Equation 3) when the initial values are in +3.

Theorem 5. The solution of model (Equation 3) is always bounded in +3 for the initial condition in +3.

Proof. Since we work the population model, it is natural that the population must be bounded due to the limitation of their biological resources, which is also known as environmental carrying capacity. Thus, the boundedness of the solution of the model (Equation 3) is also important to learn and prove. From Theorem 4, we can define a positive function as follows:

N(t)=x+myn+mzn.    (19)

For any γ > 0, the following fractional-order differential equation holds.

       CDtαN(t)+γN(t)= CDtαx(t)+mnCDtαy(t)+mnCDtαy(t)+γN(t)= (rx(1-xK)-mxz-hxc+x)      +mn(nxz-βy-δ1y-ωy2)      +mn(βy-δ2z)+γx+γmyn+γmzn= rx-rx2K-hxc+x-δ1myn-ωmy2n      -δ2mzn+γx+γmyn+γmzn rx-rx2K-δ1myn-δ2mzn+γx+γmyn+γmzn= (r+γ)x-rx2K+(γ-δ1)myn+(γ-δ2)mzn

By choosing γ < min{δ1, δ2}, we obtain

 CDtαN(t)+γN(t) (r+γ)x-rx2K=-rK(x2-(r+γ)Kxr)=-rK(x2-(r+γ)Kxr+(r+γ)2K24r2-(r+γ)2K24r2)=-rK[(x-(r+γ)K2r)2-(r+γ)2K24r2]=-rK(x-(r+γ)K2r)2+(r+γ)2K4r (r+γ)2K4r

According to Lemma 3 in Panigoro et al. [34], we apply the comparison principle and obtain

N(t)(N(0)-(r+γ)2K4γr)Eα[-γtα]+(r+γ)2K4γr.    (20)

For t → ∞, we achieve N(t)(r+γ)2K4γr, which means all solutions of model (Equation 3) with non-negative initial conditions are confined to Ω^ where

Ω^:={(x,y,z)+3 : N(t)=x(t)+my(t)n+mz(t)nσ, σ=(r+γ)2K4γr+ε, ε>0}.    (21)

4. Global dynamics

In this section, the global dynamics of model (Equation 3) are investigated. Note that all biological equilibrium points, their existence conditions, and their local stability are successfully described in Panigoro et al. [29], which can be rewritten by the following theorem.

Theorem 6. (i) The origin point Eo=(0,0,0) always exists. It is locally asymptotically stable (LAS) if r<hc.

(ii) The axial point EA=(x^,0,0) where x^ is the positive root of x2+(c-K)x+(hr-c)K=0, which has

(a) an equilibrium point if c>hr.

(b) a pair of equilibrium points if c<min{K,hr}.

Moreover, it is LAS if h<(c+x^)2rK and x^<(β+δ1)δ2βn.

(iii) The interior point EI=(x*,y*,z*) exists, if ai, i = 2, 3 in Panigoro et al. [29] satisfies the following statements.

(a) An equilibrium point in +3 if a3 < 0.

(b) Two equilibrium points in +3 if a2 < 0 and a3 > 0.

The LAS condition of EI can be seen in Theorem 4 in Panigoro et al. [29].

Note that all equilibrium points may attain local asymptotic stability with several biological conditions. Now, we will identify the biological properties to obtain globally asymptotically stable (GAS) for each equilibrium point. The analytical results are given by the following three theorems.

Theorem 7. The origin point Eo=(0,0,0) is GAS if rhc+σ.

Proof. We define the positive definite Lyapunov function as follows:

V1(x,y,z)=x+myn+mzn    (22)

By applying Lemma 3.1 in Vargas-De-León [35], we compute the α−order derivative of V1(x,y,z) along the solution of the model (Equation 3) as follows:

 CDtαV1(x,y,z) CDtαx+mnCDtαy+mnCDtαz                            =(rx(1-xK)-mxz-hxc+x)                                  +mn(nxz-βy-δ1y-ωy2)                                  +mn(βy-δ2z)                            =rx-rx2K-mxz-hxc+x+mxz                                  -βmyn-δ1myn-ωmy2n+βmyn-δ2mzn                            =rx-rx2K-hxc+x-δ1myn-ωmy2n-δ2mzn                             rx-hxc+x-δ1myn-δ2mzn.

From Equation (21), we have x ≤ σ and hence

 CDtαV1(x,y,z) rx-hxc+σ-δ1myn-δ2mzn                            =-(hc+σ-r)x-δ1myn-δ2mzn

Therefore, CDtαV1(x,y,z)0 for all (x,y,z)+3, if rhc+σ. We also find that CDtαV1(x,y,z)=0, if (x, y, z) = (0, 0, 0). This conveys that {Eo} is the only invariant set on which CDtαV1(x,y,z)=0. Obeying Lemma 4.6 in Huo et al. [20], rhc+σ obviously becomes the biological condition of Eo to reach GAS.

Theorem 8. The axial point EA=(x^,0,0) is GAS if hKcr<x^<δ2n.

Proof. We construct a positive definite Lyapunov function based on the Volterra equation as follows:

V2(x,y,z)=(x-x^-x^lnxx^)+myn+mzn.    (23)

The α-order derivative of V2(x,y,z) along the solution of the model (Equation 3) given by Lemma 3.1 in Vargas-De-León [35] is given by

 CDtαV2(x,y,z) (x-x^x)CDtαx+mCDtαyn+mCDtαzn                           =(x-x^x)(rx(1-xK)-mxz-hxc+x)                           Z+mn(nxz-βy-δ1y-ωy2)+                                 mn(βy-δ2z)                           =(x-x^)(r-rxK-mz-hc+x)+mxz                                 -mδ1yn-mωy2n-mδ2zn                           =(x-x^)(-rK(x-x^)+h(x-x^)(c+x)(c+x^)-mz)                                 +mxz-mδ1yn-mωy2n-mδ2zn                           =-rK(x-x^)2+h(x-x^)2(c+x)(c+x^)+mx^z                                 -mδ1yn-mωy2n-mδ2zn                            -rK(x-x^)2+h(x-x^)2cx^                                 +mx^z-mδ1yn-mδ2zn                           =-(rK-hcx^)(x-x^)2-mδ1yn                                 -(δ2n-x^)mz

Since hKcr<x^<δ2n, we have CDtαV2(x,y,z)0 for all (x,y,z)+3. It is also clear that CDtαV2(x,y,z)=0 if (x,y,z)=(x^,0,0). This confirms that {EA} is the only invariant set on which CDtαV2(x,y,z)=0. Therefore, EA is GAS due to Lemma 4.6 in Huo et al. [20]. This confirms the justifiability of Theorem 8.

Theorem 9. Let ΩX:={(x,y,z) : z*z<(1-mx*)my*-nσ(1+σm)my*} and h<c2rK. The interior point EI=(x*,y*,z*) is GAS in ΩX.

Proof. Consider a positive definite Lyapunov function as follows:

V3(x,y,z)=(x-x*-x*lnxx*)+mn(y-y*-y*lnxy*)+1δ2(z-z*)22z*.    (24)

By applying Lemma 3.1 in Vargas-De-León [35] and Lemma 1 in Aguila-Camacho et al. [36], we obtain the α-order derivative of V3(x,y,z) as follows:

 CDtαV3(x,y,z) (x-x*x)CDtαx+mn(y-y*y)                                  CDtαy+1δ2(z-z*z*)CDtαz                           =(x-x*x)(rx(1-xK)-mxz-hxc+x)                                 +mn(y-y*y)(nxz-βy-δ1y-ωy2)                                 +1δ2(z-z*z*)(βy-δ2z)                            =(x-x*)(-rK(x-x*)-m(z-z*)                                 +h(x-x*)(c+x)(c+x*))                                 +(y-y*)(mxzy-mx*z*y*-mω(y-y*)n)                                 +(z-z*)(yy*-1-z-z*z*)                           =-rK(x-x*)2+mz*x+mx*z                                 +h(x-x*)2(c+x)(c+x*)                                 *-my*xzy-mx*z*yy*                                 -mω(y-y*)2n+yzy*-z*yy*                                 -z+z*-(z-z*)2z*.

Applying Equation (21), we have

 CDtαV3(x,y,z) -(rK-hc2)(x-x*)2                             -mω(y-y*)2n-(z-z*)2z*                             -(1-mx*-σnmy*)z+(1+σm)z*.

Since z*z<(1-mx*)my*-nσ(1+σm)my*, we achieve

 CDtαV2(x,y,z)-(rK-hc2)(x-x*)2                           -mω(y-y*)2n-(z-z*)2z*.

Thus, CDtαV3(x,y,z)0 for all (x,y,z)+3, when h<c2rK. We also confirm that CDtαV3(x,y,z)=0 if (x, y, z) = (x*, y*, z*) and hence {EI} is the only invariant set on which CDtαV3(x,y,z)=0. Based on Lemma 4.6 in Huo et al. [20], the interior point EI is GAS in ΩX. This ends the proof.

5. Bifurcation analysis

In this section, we will study the occurrence of several phenomena namely transcritical, saddle-node, backward, and Hopf bifurcations. Two parameters are chosen, namely the harvesting rate (h) and the order of the derivative (α), as the memory index. For the analytical purpose, we define the following parameter.

h1*= cr,h2*= (c+K)2r4K,α*= 2πarctan|ζ2ζ1|.

Next, the following theorem is given for describing the occurrence of transcritical bifurcation driven by the harvesting rate (h) as the bifurcation parameter.

Theorem 10. The origin point Eo and the axial point EA exchange their stability via transcritical bifurcation when h passes through h1*.

Proof. Since the axial consists of two equilibrium points, we focus on the axial point nearest to the origin point. When h=h1*, the axial point merge with the origin point E0=EA=(0,0,0) where the eigenvalues of the Jacobian matrix are: λ1 = 0, λ2 = (β + δ1), and λ3 = −δ2. We obtain arg(λ2,3) = π > απ/2 while arg(λ1) = απ/2. This means E0=EA=(0,0,0) is non-hyperbolic. When h1*<h<(c+K)2r4K, by applying Theorems 2 and 3 in Panigoro et al. [29], E0 becomes LAS while the nearest EA becomes a saddle point. For 0<h<h1*, The origin E0 becomes unstable and the nearest EA+3 becomes unstable. This condition shows the existence of transcritical bifurcation, where h becomes the bifurcation parameter while h=h1* is the bifurcation point.

Now, the existence of saddle-node bifurcation on axial will be proven by still regarding the harvesting rate (h) as the bifurcation parameter. As a result, the following theorem is proposed.

Theorem 11. Suppose that c<min{hr,K}. The axial point EA undergoes saddle-node bifurcation when h passes through the bifurcation point h2*.

Proof. According to Theorem 1 in Panigoro et al. [29], the axial point does not exist when h>h2*. When h=h2*, a unique equilibrium point EA=(K-c2,0,0) occurs in axial where its Jacobian matrix has three eigenvalues: λ1 = 0 and λ2,3=-12[β+δ1+δ2+(β+δ1-δ2)2+2βn(K-c)]. Since |arg(λ1)| = απ/2, this axial point is non-hyperbolic. When h < h*, two axial points occurs given by EAa=(x^a,0,0) and EAb=(x^b,0,0), where x^a=K-c2+(h*-h)Kr and x^b=K-c2-(h*-h)Kr. It is easy to validate that both EAa and EAb are in +3 and have different stability. As a consequence, all the given circumstances express the occurrence of saddle-node bifurcation.

Based on Theorems 10 and 11, we obtain more global bifurcation namely backward bifurcation given by the following lemma.

Lemma 1. The model (Equation 3) undergoes backward bifurcation driven by the harvesting rate (h).

Proof. From previous theorems, the axial point EAa exists and is LAS, while Eo is unstable when h<h1*. When h1*<h<h2*, Eo becomes LAS, EAa still exists and is LAS, and unstable EAb occurs. The bistability condition is held for this interval of h which means that the convergence of the solution is very sensitive to the initial condition. Finally, those two axial points merge when h=h2* and disappear when h>h2*. This completes the proof.

Finally, we will show that the memory index in this case, that is, the order of the derivative (α), affects the dynamical behaviors of the model (Equation 3) indicated by the appearance of Hopf bifurcation around the interior point EI.

Theorem 12. Suppose the characteristic equation of the Jacobian matrix evaluated at EI can be written as λ3+ξ1λ2+ξ2λ+ξ3=0, which has a pair of complex conjugate eigenvalues λ1,2 = ζ1 ± 2, where ζ1 > 0 and one real negative eigenvalue (λ3 < 0). The model (Equation 3) undergoes a Hopf bifurcation when the order of the fractional derivative α crosses out the critical value α*=2πarctan|ζ2ζ1|.

Proof. From the earlier assumptions, we have min1i3|arg(λi)|=arctan|ζ2ζ1|. Therefore, the solution of m(α*)=α*π2-min|arg(λi)|=0 is only when α*=2πarctan|ζ2ζ1|. If we check the transversal condition: dm(α)dα|α=α*=π2 which is not equal to 0, we can assure that the sign of m(α) changes when the bifurcation parameter α passing by α*. It means that the equilibrium point EI is stable when α ∈ (0, α*) and is unstable for α* < α < 1.

6. Numerical simulations

In this section, we explore the dynamical behaviors of the model (Equation 3) numerically to support analytical findings, especially the occurrence of bifurcation. The predictor–corrector scheme given by Diethelm et al. [37] is employed. All of the parameters used in these simulations are assumptions matched with the biological conditions given by the previous analysis results. This decision was taken because this work does not specifically address an ecological case involving a particular organism.

To show the occurrence of several bifurcations driven by the harvesting rate (h), we first set the parameter as follows:

r=0.1, K=5, m=0.25, c=0.5, n=0.01,   β=0.06, δ1=0.05, δ2=0.05, ω=0.1, α=0.9.    (25)

By varying the harvesting rate in the interval 0 ≤ h ≤ 0.24, the bifurcation diagram is portrayed as in Figure 1. We have three types of dynamic behaviors around the axial point. When 0h<h1*=0.05, we have unstable origin point Eo and LAS EA. The origin point losses its stability via transcritical bifurcation when h crosses h1* and the unstable axial point EA occurs simultaneously. These dynamics are maintained for interval h1*<h<h2*=0.15125. On the other hand, the stable branch of axial point EA is preserved for 0h<h2*. The LAS point and unstable point of EA merge into the non-hyperbolic point when h=h2*. The axial point finally disappeared when h passes through h2* while the sign of Eo does not change. Thus, we have saddle-node bifurcation on axial with h2* as the bifurcation point. If we observe from a more global point of view, these interesting phenomena represent the existence of backward bifurcation marked by the occurrence of bistability condition. To show these dynamical behaviors, we choose three values of harvesting rate in each interval: h = 0.02, 0.12, and 0.18 and portray them as phase portraits and time series (see Figures 24). The interesting phenomenon called bistability is portrayed in Figure 3. Two equilibrium points LAS simultaneously impact the sensitivity of the convergence of the solution to the selection of the initial value. The two closest initial values are set which converge to the different equilibrium points. One of them convergent to the origin point and the other solution convergent to the axial point. This means, two conditions may arrive, namely the extinction of all populations and the only prey existence point. Several references show that the bistability condition occurs as the consequence of saddle-node bifurcation, see Adhikary et al. [38] and several references therein.

FIGURE 1
www.frontiersin.org

Figure 1. Bifurcation diagram driven by the harvesting rate (h) of the model (Equation 3) around the axial point using the parameter values: r = 0.1, K = 5, m = 0.25, c = 0.5, n = 0.01, β = 0.06, δ1 = 0.05, δ2 = 0.05, ω = 0.1, and α = 0.9.

FIGURE 2
www.frontiersin.org

Figure 2. Phase portrait and time series of the model (Equation 3) using parameter values: r = 0.1, K = 5, m = 0.25, c = 0.5, n = 0.01, β = 0.06, δ1 = 0.05, δ2 = 0.05, ω = 0.1, α = 0.9, and h = 0.02.

FIGURE 3
www.frontiersin.org

Figure 3. Phase portrait and time series of the model (Equation 3) using parameter values: r = 0.1, K = 5, m = 0.25, c = 0.5, n = 0.01, β = 0.06, δ1 = 0.05, δ2 = 0.05, ω = 0.1, α = 0.9, and h = 0.12.

FIGURE 4
www.frontiersin.org

Figure 4. Phase portrait and time series of the model (Equation 3) using parameter values: r = 0.1, K = 5, m = 0.25, c = 0.5, n = 0.01, β = 0.06, δ1 = 0.05, δ2 = 0.05, ω = 0.1, α = 0.9, and h = 0.18.

From the biological point of view, these numerical bifurcations describe the possibility for the prey to extinct or survive due to the change in the harvesting rate while the predator both mature and immature is always extinct (Figure 1). Three feasible conditions may happen. First, for any sufficiently small harvesting rate (0h<h1*=0.05), the prey population may maintain its existence in this ecosystem (Figure 2). Second, if the harvesting rate increases (h1*<h<h2*), two possible conditions may occur namely the extinction of prey or the viability of prey. These circumstances depend on the initial condition, where if the initial condition is quite close to the origin point, the prey will be extinct, and for the initial condition is far enough from the origin point, the prey can survive extinction (Figure 3). Third, if the harvesting rate is again increased (h>h2*), the population of prey will finally extinct and there are no population again in this ecosystem (Figure 4).

The next circumstance occurs in the interior point of the model (Equation 3), which demonstrates the influence of the order of the derivative as the memory index on the dynamical behaviors around the interior point. We set the parameter as follows:

r=0.8, K=5, m=0.25, h=0.01, c=0.08, n=0.2,   β=0.4, δ1=0.01, δ2=0.01, δ2=0.01, ω=0.1.    (26)

To identify the dynamical behaviors, we vary the values of α in the interval 0.76 ≤ α ≤ 1. As a result, we obtain the bifurcation diagram given in Figure 5. For α < α* ≈ 0.86, the interior point EI is LAS. To show this condition, we give the phase portraits by selecting α = 0.81 and α = 0.84 as given in Figures 6A, B. Nearby solution oscillates and convergent to EI. When α crosses α* ≈ 0.86, EI losses its stability via Hopf bifurcation which is indicated by the occurrence of periodic signal namely limit-cycle. The nearby solution stays away from EI and convergent to the limit-cycle. The evolution of the limit-cycle given in Figure 5 also shows that the diameter of the limit-cycle increases when alpha increases. We portray the phase portraits in Figures 6C, D to show the dynamics of solutions around EI for α = 0.87 and α = 0.9. It is shown that the densities of all populations are oscillated and finally converge to the limit cycle. The physical interpretations of Hopf bifurcation driven by α are closely related to the influence of the memory on the change of behaviors around the interior point. The stronger the influence of memory, the higher the ability of prey and predators to maintain their existence (α < α*). For less memory effect which is indicated by α > α*, all populations lose the ability to stabilize their number of individuals. The population density fluctuates according to seasonal patterns which indicates by the appearance of a limit cycle (Figures 6C, D), and the peak of each population also increases for less memory effect (Figure 7). Although the density for each population cannot tend to a constant number, in this case, the memory effect cannot cause the extinction of every population.

FIGURE 5
www.frontiersin.org

Figure 5. Bifurcation diagram driven by the order of the derivative (α) of model (Equation 3) around the axial point EI using parameter values: r = 0.8, K = 5, m = 0.25, h = 0.01, c = 0.08, n = 0.2, β = 0.4, δ1 = 0.01, δ2 = 0.01, δ2 = 0.01, and ω = 0.1.

FIGURE 6
www.frontiersin.org

Figure 6. Phase portrait of the model (Equation 3) around interior point EI using parameter values from Equation (26). (A) α = 0.81, (B) α = 0.84, (C) α = 0.87, and (D) α = 0.9.

FIGURE 7
www.frontiersin.org

Figure 7. Phase portrait of the model (Equation 3) around interior point EI using parameter values from Equation (26).

7. Conclusion

The dynamics of a predator–prey model incorporating four biological conditions, namely age structure, intraspecific competition, Michaelis–Menten type harvesting, and memory effect, have been studied. All biological validities have been presented such as the existence, uniqueness, non-negativity, and boundedness of the solution. The dynamics of the model have been explored by showing the global stability conditions for three equilibrium points namely the origin, the axial, and the interior points. We have shown that three bifurcations phenomena driven by the harvesting rate occur around the axial point namely transcritical, saddle-node, and backward bifurcations. The bistability condition exists as the impact of saddle-node bifurcation which states that the existence of prey depends on the initial condition. A bifurcation driven by the memory effect also has been shown around the interior point which is called Hopf bifurcation. All the bifurcation phenomena having an impact on the densities of the population not only may reduce their densities but also threaten the existence of several populations.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author/s.

Author contributions

All authors listed have made a substantial, direct, and intellectual contribution to the work and approved it for publication.

Funding

This research was funded by LPPM-UNG via PNBP-Universitas Negeri Gorontalo according to DIPA-UNG No. 023.17.2.677521/2021, under contract No. B/176/UN47.DI/PT.01.03/2022.

Acknowledgments

The authors intended to express our gratitude and appreciation to the editors, reviewers, and all those who have supported us in improving this article.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Deng H, Chen F, Zhu Z, Li Z. Dynamic behaviors of Lotka-Volterra predator-prey model incorporating predator cannibalism. Adv Diff Equat. (2019) 2019:359. doi: 10.1186/s13662-019-2289-8

CrossRef Full Text | Google Scholar

2. Huang Q, Lin Y, Zhong Q, Ma F, Zhang Y. The impact of microplastic particles on population dynamics of predator and prey: implication of the Lotka-Volterra model. Sci Rep. (2020) 10:1–10. doi: 10.1038/s41598-020-61414-3

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Tahara T, Gavina MKA, Kawano T, Tubay JM, Rabajante JF, Ito H, et al. Asymptotic stability of a modified Lotka-Volterra model with small immigrations. Sci Rep. (2018) 8:1–7. doi: 10.1038/s41598-018-25436-2

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Zeng X, Liu L, Xie W. Existence and uniqueness of the positive steady state solution for a Lotka-Volterra predator-prey model with a crowding term. Acta Math Sci. (2020) 40:1961–80. doi: 10.1007/s10473-020-0622-7

CrossRef Full Text | Google Scholar

5. Gui Z, Ge W. The effect of harvesting on a predator-prey system with stage structure. Ecol Modell. (2005) 187:329–40. doi: 10.1016/j.ecolmodel.2005.01.052

CrossRef Full Text | Google Scholar

6. Magnússon KG. Destabilizing effect of cannibalism on a structured predator-prey system. Math Biosci. (1999) 155:61–75. doi: 10.1016/S0025-5564(98)10051-2

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Liu C, Zhang Q, Huang J. Stability analysis of a harvested prey-predator model with stage structure and maturation delay. Math Problems Eng. (2013) 2013:329592. doi: 10.1155/2013/329592

CrossRef Full Text | Google Scholar

8. Dubey B, Agarwal S, Kumar A. Optimal harvesting policy of a prey-predator model with Crowley-Martin-type functional response and stage structure in the predator. Nonlinear Anal Model Control. (2018) 23:493–514. doi: 10.15388/NA.2018.4.3

CrossRef Full Text | Google Scholar

9. Xiao Z, Li Z, Zhu Z, Chen F. Hopf bifurcation and stability in a Beddington-DeAngelis predator-prey model with stage structure for predator and time delay incorporating prey refuge. Open Math. (2019) 17:141–59. doi: 10.1515/math-2019-0014

CrossRef Full Text | Google Scholar

10. Zhang X, Liu Z. Periodic oscillations in age-structured ratio-dependent predator-prey model with Michaelis-Menten type functional response. Physica D. (2019) 389:51–63. doi: 10.1016/j.physd.2018.10.002

CrossRef Full Text | Google Scholar

11. Zhang F, Chen Y, Li J. Dynamical analysis of a stage-structured predator-prey model with cannibalism. Math Biosci. (2019) 307:33–41. doi: 10.1016/j.mbs.2018.11.004

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Lu W, Xia Y. Periodic solution of a stage-structured predator-prey model with Crowley-Martin type functional response. AIMS Math. (2022) 7:8162–75. doi: 10.3934/math.2022454

CrossRef Full Text | Google Scholar

13. Li J, Liu X, Wei C. The impact of role reversal on the dynamics of predator-prey model with stage structure. Appl Math Model. (2022) 104:339–57. doi: 10.1016/j.apm.2021.11.029

CrossRef Full Text | Google Scholar

14. Wang W, Chen L. A predator-prey system with stage-structure for predator. Comput Math Appl. (1997) 33:83–91. doi: 10.1016/S0898-1221(97)00056-4

CrossRef Full Text | Google Scholar

15. Chow C, Hoti M, Li C, Lan K. Local stability analysis on Lotka-Volterra predator-prey models with prey refuge and harvesting. Math Methods Appl Sci. (2018) 41:7711–32. doi: 10.1002/mma.5234

CrossRef Full Text | Google Scholar

16. Xiao M, Cao J. Hopf bifurcation and non-hyperbolic equilibrium in a ratio-dependent predator-prey model with linear harvesting rate: analysis and computation. Math Comput Model. (2009) 50:360–79. doi: 10.1016/j.mcm.2009.04.018

CrossRef Full Text | Google Scholar

17. Zhang Z, Upadhyay RK, Datta J. Bifurcation analysis of a modified Leslie-Gower model with Holling type-IV functional response and nonlinear prey harvesting. Adv Diff Equat. (2018) 2018:127. doi: 10.1186/s13662-018-1581-3

CrossRef Full Text | Google Scholar

18. Mulheisen M, Allen C, Parr CS. Lycaon Pictus. (On-line), Animal Diversity Web. (2002). Available online at: https://animaldiversity.org/accounts/Lycaon_pictus/ (accessed on November 13, 2022)

19. IUCN SSC Antelope Specialist Group. Aepyceros melampus. The IUCN Red List of Threatened Species 2016: eT550A50180828. (2016). Available online at: https://www.iucnredlist.org/species/550/50180828 (accessed on November 13, 2022)

20. Huo J, Zhao H, Zhu L. The effect of vaccines on backward bifurcation in a fractional order HIV model. Nonlinear Anal. Real World Appl. (2015) 26:289–305. doi: 10.1016/j.nonrwa.2015.05.014

CrossRef Full Text | Google Scholar

21. Moustafa M, Mohd MH, Ismail AI, Abdullah FA. Stage structure and refuge effects in the dynamical analysis of a fractional order Rosenzweig-MacArthur prey-predator model. Progr Fract Different Appl. (2019) 5:49–64. doi: 10.18576/pfda/050106

CrossRef Full Text | Google Scholar

22. Rahmi E, Darti I, Suryanto A, Trisilowati, Panigoro HS. Stability analysis of a fractional-order leslie-gower model with allee effect in predator. J Phys. (2021) 1821:012051. doi: 10.1088/1742-6596/1821/1/012051

CrossRef Full Text | Google Scholar

23. Owolabi KM. Dynamical behaviour of fractional-order predator-prey system of Holling-type. Discrete Contin Dyn Syst S. (2020) 13:823–34. doi: 10.3934/dcdss.2020047

CrossRef Full Text | Google Scholar

24. Barman D, Roy J, Alam S. Modelling hiding behaviour in a predator-prey system by both integer order and fractional order derivatives. Ecol Inform. (2022) 67:101483. doi: 10.1016/j.ecoinf.2021.101483

CrossRef Full Text | Google Scholar

25. Ghanbari B, Djilali S. Mathematical and numerical analysis of a three-species predator-prey model with herd behavior and time fractional-order derivative. Math Methods Appl Sci. (2020) 43:1736–52. doi: 10.1002/mma.5999

CrossRef Full Text | Google Scholar

26. Yousef FB, Yousef A, Maji C. Effects of fear in a fractional-order predator-prey system with predator density-dependent prey mortality. Chaos Solitons Fractals. (2021) 145:110711. doi: 10.1016/j.chaos.2021.110711

CrossRef Full Text | Google Scholar

27. Ghosh U, Pal S, Banerjee M. Memory effect on Bazykin's prey-predator model: stability and bifurcation analysis. Chaos Solitons Fractals. (2021) 143:110531. doi: 10.1016/j.chaos.2020.110531

CrossRef Full Text | Google Scholar

28. Panigoro HS, Rahmi E, Achmad N, Mahmud SL, Resmawan R, Nuha AR. A discrete-time fractional-order rosenzweig-macarthur predator-prey model involving prey refuge. Commun Math Biol Neurosci. (2021) 2021:1–19. doi: 10.28919/cmbn/6586

CrossRef Full Text | Google Scholar

29. Panigoro HS, Resmawan R, Sidik ATR, Walangadi N, Ismail A, Husuna C. A fractional-order predator-prey model with age structure on predator and nonlinear harvesting on prey. Jambura J Math. (2022) 4:355–66. doi: 10.34312/jjom.v4i2.15220

CrossRef Full Text | Google Scholar

30. Mahata A, Paul S, Mukherjee S, Das M, Roy B. Dynamics of caputo fractional order SEIRV epidemic model with optimal control and stability analysis. Int J Appl Comput Math. (2022) 8:1–25. doi: 10.1007/s40819-021-01224-x

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Yavuz M, Sene N. Stability analysis and numerical computation of the fractional predator-prey model with the harvesting rate. Fractal Fract. (2020) 4:35. doi: 10.3390/fractalfract4030035

CrossRef Full Text | Google Scholar

32. Maji C. Impact of fear effect in a fractional-order predator-prey system incorporating constant prey refuge. Nonlinear Dyn. (2022) 107:1329–42. doi: 10.1007/s11071-021-07031-9

CrossRef Full Text | Google Scholar

33. Barman D, Roy J, Alrabaiah H, Panja P, Mondal SP, Alam S. Impact of predator incited fear and prey refuge in a fractional order prey predator model. Chaos Solitons Fractals. (2021) 142:110420. doi: 10.1016/j.chaos.2020.110420

CrossRef Full Text | Google Scholar

34. Panigoro HS, Suryanto A, Kusumahwinahyu WM, Darti I. Dynamics of a fractional-order predator-prey model with infectious diseases in prey. Commun Biomath Sci. (2019) 2:105. doi: 10.5614/cbms.2019.2.2.4

CrossRef Full Text | Google Scholar

35. Vargas-De-León C. Volterra-type Lyapunov functions for fractional-order epidemic systems. Commun Nonlinear Sci Num Simulat. (2015) 24:75–85. doi: 10.1016/j.cnsns.2014.12.013

CrossRef Full Text | Google Scholar

36. Aguila-Camacho N, Duarte-Mermoud MA, Gallegos JA. Lyapunov functions for fractional order systems. Commun Nonlinear Sci Num Simulat. (2014) 19:2951–7. doi: 10.1016/j.cnsns.2014.01.022

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Diethelm K, Ford NJ, Freed AD. A predictor-corrector approach for the numerical solution of fractional differential equations. Nonlinear Dyn. (2002) 29:3–22. doi: 10.1023/A:1016592219341

CrossRef Full Text | Google Scholar

38. Adhikary PD, Mukherjee S, Ghosh B. Bifurcations and hydra effects in Bazykin's predator-prey model. Theor Populat Biol. (2021) 140:44–53. doi: 10.1016/j.tpb.2021.05.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Keywords: bifurcation, age structure, intraspecific competition, harvesting, memory effect

Citation: Panigoro HS, Rahmi E and Resmawan R (2023) Bifurcation analysis of a predator–prey model involving age structure, intraspecific competition, Michaelis–Menten type harvesting, and memory effect. Front. Appl. Math. Stat. 8:1077831. doi: 10.3389/fams.2022.1077831

Received: 23 October 2022; Accepted: 25 November 2022;
Published: 06 January 2023.

Edited by:

Bapan Ghosh, Indian Institute of Technology Indore, India

Reviewed by:

Huda Abdul Satar, University of Baghdad, Iraq
Prabir Panja, Haldia Institute of Technology, India

Copyright © 2023 Panigoro, Rahmi and Resmawan. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Hasan S. Panigoro, yes hspanigoro@ung.ac.id

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.