^{*}

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

This article was submitted to Mathematical Biology, a section of the journal Frontiers in Applied Mathematics and Statistics

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.

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.

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. [

In Wang and Chen [_{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

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 [

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 [

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. [

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. [_{i}(_{1}, _{3}, ∥_{4}, ∥_{5}, and _{2} = β + δ_{1} + ω(_{3} + _{4}), and _{3} = δ_{2}. Therefore, we conclude that _{i}(_{i} < 1, then _{i}(

Theorem 1. The kernel _{i}(_{i} < 1,

Next, Equation (5) can be written as follows:
_{0}(_{0}(_{0}(

Theorem 2. The solution of model (Equation 3) has a solution under the condition if we have _{1} such that

_{i}, _{in}, _{i∞} → 0 ∀_{1}, we have
_{i∞}∥ → 0 ∀

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 _{1}(_{1}(_{1}(_{1}(_{1}(_{1}(_{1}(_{1}(_{1}(

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

Theorem 4. If the initial condition in

Theorem 5. The solution of model (Equation 3) is always bounded in

_{1}, δ_{2}}, we obtain

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. [

Theorem 6.

The origin point

The axial point

an equilibrium point if

a pair of equilibrium points if

Moreover, it is LAS if

The interior point _{i},

An equilibrium point in _{3} < 0.

Two equilibrium points in _{2} < 0 and _{3} > 0.

The LAS condition of

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

Theorem 8. The axial point

Theorem 9. Let _{X}.

^{*}, ^{*}, ^{*}) and hence _{X}. This ends the proof.

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 (

Theorem 10. The origin point

_{1} = 0, λ_{2} = (β + δ_{1}), and λ_{3} = −δ_{2}. We obtain arg(λ_{2,3}) = π > απ/2 while arg(λ_{1}) = απ/2. This means

Now, the existence of saddle-node bifurcation on axial will be proven by still regarding the harvesting rate (

Theorem 11. Suppose that

_{1} = 0 and _{1})| = απ/2, this axial point is non-hyperbolic. When ^{*}, two axial points occurs given by

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 (

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

Theorem 12. Suppose the characteristic equation of the Jacobian matrix evaluated at _{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

^{*}. It means that the equilibrium point ^{*}) and is unstable for α^{*} < α < 1.

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. [

To show the occurrence of several bifurcations driven by the harvesting rate (

Bifurcation diagram driven by the harvesting rate (_{1} = 0.05, δ_{2} = 0.05, ω = 0.1, and α = 0.9.

Phase portrait and time series of the model (Equation 3) using parameter values: _{1} = 0.05, δ_{2} = 0.05, ω = 0.1, α = 0.9, and

Phase portrait and time series of the model (Equation 3) using parameter values: _{1} = 0.05, δ_{2} = 0.05, ω = 0.1, α = 0.9, and

Phase portrait and time series of the model (Equation 3) using parameter values: _{1} = 0.05, δ_{2} = 0.05, ω = 0.1, α = 0.9, and

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 (

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:
^{*} ≈ 0.86, the interior point ^{*} ≈ 0.86, ^{*}). 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 (

Bifurcation diagram driven by the order of the derivative (α) of model (Equation 3) around the axial point _{1} = 0.01, δ_{2} = 0.01, δ_{2} = 0.01, and ω = 0.1.

Phase portrait of the model (Equation 3) around interior point

Phase portrait of the model (Equation 3) around interior point

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.

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

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

This research was funded by LPPM-UNG

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

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.

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.