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 [6–13].

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: (1)dxdt= rx(1-xK)-mxz,dydt= nxz-βy-δ1y,dzdt= βy-δ2z. 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, δ₁ and δ₂ 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 (2)dxdt= rx(1-xK)-mxz-hxc+x,dydt= nxz-βy-δ1y-ωy2,dzdt= βy-δ2z. 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: (3)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). 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: (4) CDtαf(t)=1Γ(1-α)∫0tf′(τ)(t-τ)αdτ, 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.

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. (5)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τ, Now, we will show that the kernels F_i(x, y, z), i = 1, 2, 3, satisfy the Lipschitz condition. For ∥.∥ is the Euclidean norm, we suppose that ∥x(t)∥ ≤ a₁, ∥x̄(t)∥≤a2, ∥y(t)∥ ≤ a₃, ∥ȳ(t)∥ ≤ a₄, ∥z(t)∥ ≤ a₅, 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̄))||≤ r∥x-x̄∥+(a1+a2)rK∥x-x̄∥+a5m∥x-x̄∥+hc∥x-x̄∥= g1∥x-x̄∥,      ∥F2(x,y,z)-F2(x,ȳ,z)∥= ∥[nxz-βy-δ1y-ωy2]-[nxz-βȳ-δ1ȳ-ωȳ2]∥= ∥-β(y-ȳ)-δ1(y-ȳ)-ω(y+ȳ)(y-ȳ)∥≤ β∥y-ȳ∥+δ1∥y-ȳ∥+ω(a3+a4)∥y-ȳ∥= g2∥y-ȳ∥,      ∥F3(x,y,z)-F3(x,y,z̄)∥= ∥[βy-δ2z]-[βy-δ2z̄]∥≤ g3∥z-z̄∥, where g1=r+(a1+a2)rK+a5m+hc, g₂ = β + δ₁ + ω(a₃ + a₄), and g₃ = δ₂. Therefore, we conclude that F_i(x, y, z), i = 1, 2, 3, satisfy the Lipschitz condition. Furthermore, it is clear that if 0 ≤ g_i < 1, then F_i(x, y, z) are contractions for i = 1, 2, 3. Therefore, the following theorem is obtained.

Theorem 1. The kernel F_i(x, y, z), i = 1, 2, 3 satisfy the Lipschitz condition and contractions if 0 < g_i < 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 x₀(t) = x(0), y₀(t) = y(0), and z₀(t) = z(0). Therefore, we have (6)φ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τ, 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 (7)∥φ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τ||. From Theorem 1, we have that the kernel satisfy the Lipschitz condition and hence Equation (7) becomes ∥xn(t)-xn-1(t)∥≤g1Γ(α)∫0t∥xn-1-xn-2∥(t-τ)α-1 dτ,∥yn(t)-yn-1(t)∥≤g2Γ(α)∫0t∥yn-1-yn-2∥(t-τ)α-1 dτ,∥zn(t)-zn-1(t)∥≤g3Γ(α)∫0t∥zn-1-zn-2∥(t-τ)α-1 dτ. The last inequality gives (8)∥φ1n(t)∥≤g1Γ(α)∫0t∥φ1n-1(τ)∥ dτ,∥φ2n(t)∥≤g2Γ(α)∫0t∥φ2n-1(τ)∥ dτ,∥φ3n(t)∥≤g3Γ(α)∫0t∥φ3n-1(τ)∥ dτ. 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 t₁ 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 F_i, i = 1, 2, 3 satisfy the Lipschitz condition. According to Equation (8), we obtain (9)∥φ1n(t)∥≤ ∥x0(t)∥[t1g1Γ(α+1)]n,∥φ2n(t)∥≤ ∥y0(t)∥[t1g2Γ(α+1)]n,∥φ3n(t)∥≤ ∥z0(t)∥[t1g3Γ(α+1)]n, 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 (10)x(t)-x(0)=xn(t)-Q1n(t),y(t)-y(0)=yn(t)-Q2n(t),z(t)-z(0)=zn(t)-Q3n(t). where Q_in, i = 1, 2, 3 are the remaining terms. Furthermore, the given terms would be demonstrated to hold Q_i∞ → 0 ∀i = 1, 2, 3. Denoting that (11)∥Q1n∥(t)≤ ∥1Γ(α)∫0t(F1(x,y,z)-F1(xn-1,y,z)) dτ∥                ≤ 1Γ(α)∫0t∥F1(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Γ(α)∫0t∥F2(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Γ(α)∫0t∥F3(x,y,z)-F3(x,y,zn-1)∥ dτ                ≤ tg3Γ(α+1)∥z-zn-1∥. By applying a recursive pattern, we acquire (12)∥Q1n(t)∥≤[tΓ(α+1)]n+1g1nk,∥Q2n(t)∥≤[tΓ(α+1)]n+1g2nk,∥Q3n(t)∥≤[tΓ(α+1)]n+1g3nk. At the point t₁, we have (13)∥Q1n(t)∥≤[t1Γ(α+1)]n+1g1nk,∥Q2n(t)∥≤[t1Γ(α+1)]n+1g2nk,∥Q3n(t)∥≤[t1Γ(α+1)]n+1g3nk. By considering the results of Theorem 1 and applying n → ∞ to both sides, we have ∥Q_i∞∥ → 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 x₁(t), y₁(t), and z₁(t). Then, we have (14)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τ. Applying the norm on both sides, we achieve (15)∥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τ. 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. (16)∥x(t)-x1(t)∥(1-tg1Γ(α+1))≤ 0,∥y(t)-y1(t)∥(1-tg2Γ(α+1))≤ 0,∥z(t)-z1(t)∥(1-tg3Γ(α+1))≤ 0. As a result, we achieve ∥x(t) − x₁(t)∥ = 0, ∥y(t) − y₁(t)∥ = 0, and ∥z(t) − z₁(t)∥ = 0, which impact x(t) = x₁(t), y(t) = y₁(t), and z(t) = z₁(t). Then, we finally give the following theorem.

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

In this section, we will show that for any initial condition is in ℝ+3 where ℝ+3:={(x,y,z) : x≥0, y≥0, z≥0, 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 (17){x(t)    >0, when 0≤t<t^,x(t^)    =0, when t=t^x(t^+)<0, when t=t^+. From the first equation in Equation (3) along with Equation (17), we have (18) CDtαx(t)|x(t^)=0. 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: (19)N(t)=x+myn+mzn. 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{δ₁, δ₂}, 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 (20)N(t)≤(N(0)-(r+γ)2K4γr)Eα[-γtα]+(r+γ)2K4γr. 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 (21)Ω^:={(x,y,z)∈ℝ+3 : N(t)=x(t)+my(t)n+mz(t)n≤σ, σ=(r+γ)2K4γr+ε, ε>0}.

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.

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

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 r≤hc+σ.

Proof. We define the positive definite Lyapunov function as follows: (22)V1(x,y,z)=x+myn+mzn 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 r≤hc+σ. 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], r≤hc+σ 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: (23)V2(x,y,z)=(x-x^-x^lnxx^)+myn+mzn. 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.