A dynamic analysis of a prey–predator population model with a nonlinear harvesting rate

Nadia Mohdeb (Applied Mathematics Laboratory, University of Bejaia, Bejaia, Algeria)

Arab Journal of Mathematical Sciences

ISSN: 1319-5166

Article publication date: 28 February 2023




In this article, the author discusses dynamical behaviors of a prey-predator population model with nonlinear harvesting rate and offers a mathematical analysis of the model.


The design is by using modelization of populations interaction, qualitative theory of ordinary différential equations, bifurcations analysis, invariant center manifolds theory and Dulac's criterion.


The author studies the stability of solutions and the existence of periodic solutions in the model, and proves the existence of some invariant sets and the production of a transcritical together with a saddle-node bifurcation.

Practical implications

The author studies the effects of harvesting on the persistence and extinction properties and its influence in the perspectives of economic views.


The authors considers a predator–prey model with a new nonlinear form of harvesting rate. The author’s intention is to make conceptual adjustments to a well-known predator–prey model in order to incorporate the effects of harvesting.



Mohdeb, N. (2023), "A dynamic analysis of a prey–predator population model with a nonlinear harvesting rate", Arab Journal of Mathematical Sciences, Vol. ahead-of-print No. ahead-of-print. https://doi.org/10.1108/AJMS-03-2022-0052



Emerald Publishing Limited

Copyright © 2023, Nadia Mohdeb


Published in the Arab Journal of Mathematical Sciences. Published by Emerald Publishing Limited. This article is published under the Creative Commons Attribution (CC BY 4.0) license. Anyone may reproduce, distribute, translate and create derivative works of this article (for both commercial and non-commercial purposes), subject to full attribution to the original publication and authors. The full terms of this license may be seen at http://creativecommons.org/licences/by/4.0/legalcode

1. Introduction

Many researchers are interested to the dynamic of predator-prey interactions models, and have investigated the processes that affect it [1–14]. The interaction between the predator and prey species can be modeled by the classical Lotka–Volterra model [15]. When the prey species obeys the logistic growth rate, this model becomes

where x represents the size of prey population and y represents the size of predator population. Parameters a, r, k, c and d are assumed to be positive values, where

a presents the rate of predation, such that in the presence of the predator, the prey species decreases at a rate proportional to the functional response ax;

c denotes the factor of the efficiency of predation which divides a maximum per capita birth rate of the predators into a maximum per capita consumption rate;

d is the death rate of the predator which decreases exponentially, in the absence of the preys;

r is the maximum specific growth rate of the prey which grows logistically in the absence of the predator species and

k denotes the environmental carrying capacity with which the prey grows logistically in the absence of the predation.

The persistence and extinction properties of a large marine ecosystem are not only dependent on the interspecific interaction between the marine species but it is also dependent on the external environment factors. To enrich model (1), many researchers modify the nonlinear functional response function, adding some other elements like pollution, toxicity, harvesting, age of the species, refuge, etc. [2, 5–10].

In the recent decades, human ambitions seem to be endless, and the global modernization and urbanization foment the increasing human demands in searching more food and resources. However, the ongoing civilization and urbanization have expedited the process of exploitation of natural resources and altered the environment. Humans harvest the preys and the predators for the purpose of private consumption and commercial reasons. Furthermore, in the recent years, there is a substantial growth of global concerns on the impact of harvesting activities on organisms and ecosystems.

Several forms and types of harvesting in prey–predator models are already being studied; researchers have added terms to the prey or predator density like constants [16]; linear functions; functions that are linear if the density of the predator is bellow a switched value and constant otherwise [10]; functions of the form QiFixi/Qi+Fixi,i=1,2 where Qi and Fi are constants and xi, i = 1, 2 are the sizes of prey and predator populations, respectively [17], and recently by incorporating the Heaviside function [18]. The most common one of these harvesting forms is a nonzero constant or a linear harvesting rate. Models studied with a linear or constant harvesting rate exhibit far richer and complex dynamics compared to the models with no harvesting. Hu and Cao [6] argue that nonlinear harvesting is more realistic and reasonable than modeling constant-yield harvesting and constant effort harvesting. They consider a predator–prey system with a nonlinear Michaelis–Menten type of predator harvesting and demonstrate the dynamical complexity of the system with this type of harvesting effect.

Apparently, the above types of harvesting rates have their own advantages as well as disadvantages in the award of the harvest in the real world. This background serves as the motivation for the present paper; in this work, we consider a predator–prey model with a new nonlinear form of harvesting rate. Our intention is to make conceptual adjustments to a well-known predator–prey model in order to incorporate the effects of harvesting. In the present work, we take the assumption of a prey and a predator dependent harvesting rate one step further to incorporate biological and economic realism in our modeling framework. This kind of nonlinear harvesting is more realistic and reasonable than the model with constant-yield harvesting and constant effort harvesting. It can be thought as a supplement to existing literature on the dynamics of this system, since there is little literature involved in nonlinear type harvesting for the system up to now. Taking system (1) as our baseline model, we assume that harvesting takes place, and both the prey and the predator are under harvesting. The effect of harvesting on the species is assumed to be different and is examined to make further speculations on the persistence and extinction properties. We introduce harvesting nonlinear functions H1x of the prey and H2y of the predator to model (1) for discussing its dynamical characteristics. We then consider the following differential system

where α and β are positive parameters related to the harvesting effort. We investigate the existence and stability of multiple equilibria, bifurcations and study the effects of harvest on the dynamics of the predator–prey model (2). The model possesses a varied dynamic.

This paper is organized as follow. In section 2, we provide the formulation of the model which we will study here. Basic results are given in section 3; we present the boundedness and the invariance of the model considered. We then study the existence of multiple equilibria, their types and their local and global stability, the permanence and persistence of the model, bifurcations and existence of periodic solutions, especially existence of limit cycles for the model. In section 4, we give numerical simulations to illustrate the established results. Concluding remarks are presented in section 5.

2. Model formulation

Considering prey–predator model (2), we suppose that harvesting function of the prey in system (2) is H1x=x2 and harvesting of the predator is cubic H2y=y3. Moreover, the prey equation is an extended version of the logistic equation. Model (2) is

The harvesting function of the prey is quadratic and can be added to the logistic term, and the first equation becomes x.=x(rxγ)axyαx2, where γ is a real parameter, but because of the great biological interest of harvest functions and as the objective of this work is to study and examine the effects of the parameter alpha, related to harvesting effort, in the widely studied model (1), the logistic and harvesting terms are written separately.

We can show that solutions of system (3) starting from positive initial conditions are positive and forward bounded. Let

Then, D=x,yR×R; 0<x<x0,0<y<1βcax0d is a positive invariant set for system (3).

3. Preliminary results

3.1 Existence of equilibria

In this section, we inspect the existence of all equilibria of model (3). Their existence and number depend principally on the harvesting rates. Because of their biological meaning, we will only study and explore the positive equilibria. We denote α1=α+rk, d1=dcarα1 and Δ=c2a4α124βd1.

Theorem 3.1.

According to the values of the harvesting parameter β and the mortality rate d of the predator species in the absence of the prey, system (3) has the following coexistence positive equilibria:

In the case β<c2a44d1α12 and d1 > 0, or in the case d1 < 0, system (3) admits three positive equilibria: 0,0, rα1,0 and x2,y2 where

x2=rα1aα1y2 and y2=ca2α1+Δ2β.
If β=c2a44d1α12, system (3) has two equilibria: 0,0 and rα1,0. When β>c2a44d1α12 and d1 ≥ 0 there exist two positive equilibria: 0,0 and rα1,0.

Proof. An equilibrium point of system (3) satisfies

As d is positive, we get by equations x = 0 or y = 0, two equilibria 0,0 and rα1,0. On the other hand, equation
gives x=rα1aα1y. Solutions of equation
are given based on values of d1 and β compared to that of c2a44d1α12. More precisely, it follows from equations (5) and (6): two equilibria x1,y1 and x2,y2 either if β<c2a44d1α12 and d1 > 0 or if d1 < 0; one equilibrium point x1,y1 if d1 = 0, where
(7)x1=rα1aα1y1, y1=ca2α1Δ2β; 
one equilibrium point x1,ca22α1β if β=c2a44d1α12; and none if β>c2a44d1α12 and d1 > 0.
Remark 3.1.

The origin represents an equilibrium when both the prey and predator population die out and extinct, and the point rα1,0 represents the equilibrium when the prey population survives in the absence of the predator population. The nonlinear harvest strategy applied to model (3) leads to the predator extinct when the harvesting rate conditions satisfy βc2a44d1α12 and d1 ≥ 0. In all other cases, the prey and predator coexist.

3.2 Nature and stability of equilibria

In this section, we provide the nature of positive equilibria of system (3), and their stability will be discussed.

First, we can show that the equilibrium point rα1,0 is a saddle, locally unstable, if d1 < 0, and a locally stable node when d1 > 0. However, if d1 = 0, then rα1,0 is a degenerate equilibrium point, and further study is needed. In this case, in the following theorem we will study the stability of the two equilibria, origin and rα1,0 and describe the behavior of the solutions of (3) near this last equilibrium point.

Theorem 3.2.

Equilibrium (0,0) is unstable for system (3). If d1 = 0, equilibrium point rα1,0 is a saddle-node, with two hyperbolic sectors and a parabolic one.

There exists an invariant strong unstable manifold Wu tangent at point rα1,0 to the x-axis on which the behavior of system (3) is repulsive.

On one side of Wu, for each point x,y, there exists a center manifold Wcx,y tangent at point rα1,0 to the line y=α1ax+ra, and on the other side, all the center manifolds coincide and are tangent at point rα1,0 to the same line.

Proof. The eigenvalues of the Jacobian matrix associated with system (3) at the equilibrium point 0,0 are r > 0 and d<0, and then the origin is an unstable saddle.

Now we suppose that d1 = 0. First we translate the equilibrium point rα1,0 to the origin by the change of variables X=xrα1 and Y = y. System (3) becomes

By using the transvection u=X+rα1Y, v = Y and by reversing the time τ = −t, we get
where (′) denotes ddt. By exchanging the roles of u and v in system (8) by setting x¯=v, y¯=u, we obtain
We denote Ax,y=caxy+ca2α1x2+βx3. Equation y¯=0 gives
The solutions of equation (10) depend on the sign of
Because x¯,y¯ is supposed to be in the neighborhood of the origin, then since x¯ is smallest than r2, we obtain Δ->0 and equation (10) has two solutions:
We can show that y¯1 is excluded and y¯=y¯2. By the third-order Taylor development, we get
Since ca2α1>0, following the similar arguments in Ref. [19] which is based on removing of flat terms in order to get the C normal forms by using the homotopic method, we conclude that origin is a saddle-node for system (9). It follows from the stable manifold theorem [20, 21] that
  1. there exists an invariant strong unstable manifold Wu tangent at origin to the y¯-axis such that on Wu, system (9) is analytically conjugate to y¯=ry¯ and the behavior is repulsive;

  2. for any x¯0,y¯0, if x¯0>x¯Wu, where x¯Wu is the abscissa of Wu at the point whose ordinate is y¯0, it passes a center manifold Wc tangentially at origin to the x¯-axis. On Wc, system (9) is C-conjugate to x¯=x¯2+βc2a4rα12x¯3. If x¯0<x¯Wu, all the center manifolds coincide and are tangent at origin to x¯-axis and

  3. system (9) is C-conjugate to

and C0-conjugate to
Going back to the coordinates (x, y) by: inverting transvection; exchanging the roles of u and v; translation and by change time, we obtain the result.
Remark 3.2.

In the proof of theorem 3.2, more details are given about the behavior of solutions near the equilibrium point rα1,0.

Because that the coordinates y2 is positive if and only if d1 ≤ 0, equilibrium point x2,y2 is not biological meaningful if d1 > 0; we note that x2 remains positive for all values of the parameters. We will show in the following theorem that stability and nature of equilibrium x2,y2 depend essentially upon of the values of the mortality rate d of the predator and of the harvesting parameters α and β.

Theorem 3.3.

The equilibrium point x2,y2 is locally unstable if d1 ≥ 0, and it is either a linear center, a focus or a node if d1 < 0; it is locally stable if d1 < 0 and 1caα10.

Proof. We assume that Δ ≠ 0. The Jacobian matrix of system (3) at a point x,y, such that x=rα1aα1y, is given by

We first assume that d1 < 0. We have
We denote Ix,y=2ca2α1y+d1+3βy2. Its sign depends on that of Δ1c2a4α123βd1. We have Δ1 = Δ + βd1, it is positive. This means that Ix,y has two roots:
y¯=c2a4α12Δ23β and y¯¯=c2a4α12+Δ23β,
where Δ2c2a4α123βd1. We will compare y2 with respect to y¯¯. Suppose that y2<y¯¯. Then
By squaring, since ca2α1+Δ=2βy2, the left side of (11) is positive, and we get Δ<ca2α1Δ2. It follows by some calculus that −4βd1 < 0 which is absurd, and then Ix2,y2>0.

Moreover, we can show that ray2 > 0; otherwise, Δ<ca2α1+2rβa, and by squaring and simplifications, we get d<r2βa2<0 which contradicts the assumption. Hence x2,y2 is a node, a center or a focus.

To study the stability of equilibrium x2,y2, we look at

and it is easy to see that if 1caα10, then x2,y2 is stable.

If d1 = 0, then x2,y2=rα1,0 and is, by theorem 3.2, a saddle-node.

In the following result, we use Dulac's criterion [22] to show the nonexistence of limit cycles in system (3).

Theorem 3.4.

System (3) does not have a periodic solution and then it does not have a limit cycle in

D+=(x,y)R×R; x>0,y>0.

Proof. Consider a real-valued continuously differentiable Dulac function

for x > 0, y > 0 and denote
We have
Since x > 0 and y > 0, then (hf)xx,y+(hg)yx,y is not identically zero and is of one sign. Therefore, system (3) does not have a closed orbit in D+.

We can now state results about the global stability of equilibria of model (3).

Theorem 3.5.

1) If d1 < 0, then the positive equilibrium point x2,y2 is asymptotically stable; it is globally asymptotically stable in D+.

2) If d1 > 0, then the equilibrium point rα1,0 is globally asymptotically stable in D+R+×0.

Proof. It follows from theorem 3.4 that when d1 < 0, the equilibrium point x2,y2 cannot be a center for system (3). By using in addition theorem 3.3, we deduce that x2,y2 is locally asymptotically stable.

Moreover, we know by theorem 3.1 that there is not any other equilibrium of system (3) than x2,y2 in D+, which is positively invariant for system (3). It follows from theorem 3.4 that any trajectory of system (3) in D+ tends to x2,y2 when t tends to positive infinity, and the first part of the theorem is shown.

We know that if d1 > 0, then equilibrium point rα1,0 is a stable node. It is, except the origin which is instable, the unique equilibrium point of system (3) in D+R+×0, which is positively invariant for system (3); the direction of vector fields allows us to conclude the second part of the theorem.

Theorem 3.6.

System (3) is permanent and persistent if and only if d1 < 0.

Proof. We consider in the case d1 < 0, the points:

  • Ax0,y0 a point of the curve y=1βcaxd such that dca<x0<x2;

  • B the point of intersection of the horizontal line passing through A with the line y=(rakαa)x+ra;

  • E the point of intersection of the vertical line through B with the curve y=1βcaxd;

  • M the point of intersection of the horizontal line through E with the line y=(rakαa)x+ra;

  • N the point of intersection of the vertical line through M with the curve y=1βcaxd and

  • C the part of the curve y=1βcaxd connecting the points M and A.

Let D1 be the region of the plan delimited by the straight segments AB, BE, EM and MN and by the curved line C. We show that D1 is positively invariant compact for system (3). Since in the case d1 < 0, the point x2,y2 is the unique equilibrium in the first quadrant of the plan and using, according to theorem 3.5, the fact that x2,y2 is globally asymptotically stable we can then conclude.

3.3 Bifurcation results and analysis

We present in this part the different bifurcation diagrams of multiple equilibria of system (3). We examine the effects of the harvesting parameter on the prey and predator species by valuing the parameters. We will theoretically choose the parameters used in the simulations according to the criterion: check appropriately the obtained analytical results in each case and show numerically the acuteness of the obtained stability conditions. Note that model (3) exhibits a forward bifurcation.

Example 3.1.

Choosing d = 0.2, r = 0.1, k = 0.5, c = 1.8, a = 0.5 and β = 0.1, we get d1 ≤ 0 for α ≤ 0.25. The two equilibria (x2, y2) and (rα1,0) coincide at the forward bifurcation value α = 0.25, as shown in Figure 1.

Example 3.2.

Using the following parameter values d = 0.2, r = 0.1, k = 0.5, c = 1.8, a = 0.5 and β = 0.1, we get d1 ≥ 0 and β>c2a44d1α12 for 0.25 ≤ α (Figure 2). The positive equilibrium does not exist because y − 2 becomes negative (x2 is still positive).

If α is less than the value giving d1 = 0, an unstable equilibrium and a stable one exist in model (3) and the prey and the predator coexist (Figure 1). The density of prey and predator species decrease severely in result of being harvested showing that the harvesting activities have a great impact on the populations. If the harvesting rate parameter α exceeds this value, system (3) does not have positive equilibria (Figures 1 and 2), and there is a predator extinction. The predator population encounter extinction due to the high level of harvesting activities on the predator species, and only the prey population can persist; the prey population increases, this is because the increasing of harvesting activities on predator suppress the predator population and consequently reduces the predation activities on the prey population.

The model exhibits a transcritical bifurcation. This bifurcation is ecologically important and will lead to the potentially dramatic variation of the system dynamics. Existence of such kinds of bifurcations indicates that the over exploitation of resource will cause extinction of the species. The predators tend to extinction for some values of harvesting rate while the predators and prey will coexist for some other values of harvesting rate, which is one of the most exciting features in the ecosystem.

4. Numerical simulation

We give in this section some interesting numerical examples for system (3); to illustrate better the behavior of the trajectories near the equilibria, we represent the phase portraits in the half-plane (x,y)R×R; x>0. For the graphical representations, we consider different contexts by considering the different values of the parameters of the model.

Example 4.1.

Setting parameter values d = 0.2, r = 0.1, k = 0.5, c = 1.8, a = 0.5, β = 0.1 and α = 0.1, we get d1 = −0.1 < 0. The equilibrium point (x2, y2) exists and is positive and globally asymptotically stable (Figure 3). In this case, (rα1,0) is an unstable equilibrium point.

Example 4.2.

Setting d = 0.2, r = 0.1, k = 0.5, c = 1.8, a = 0.5, β = 0.1 and α = 0.25, we get d1 = 0. The equilibrium point (x2, y2) coincides with the equilibrium (rα1,0) giving a saddle-node point and a forward (a transcritical) bifurcation occurs (Figure 4); in this case, there is no other positive equilibria. A forward bifurcation diagram is illustrated in Figures 1 and 2. The center manifold of the saddle-node point divides the plan into two regions. The region above is the basin of attraction for (rα1,0).

Example 4.3.

Using the following parameter values d = 0.2, r = 0.1, k = 0.5, c = 1.8, a = 0.5, β = 0.1 and α = 0.7, we obtain d1 = 0.1 > 0. The positive equilibrium point (x2, y2) does not exist anymore, and there is no positive equilibria other than (rα1,0); this last equilibrium is stable. The phase portrait of model (3) is shown in Figure 5.

5. Conclusion

In this model, we discussed the effects of nonlinear harvesting in a predator–prey system in which both the species are harvested. Qualitative analysis reveals that harvesting plays an important role in determining the dynamics and bifurcations of the model. The equilibria of model (3) are examined and the stability is discussed, and the parameters α and β in nonlinear harvesting term affect the number and stability of equilibria. We showed that there occurs an interesting bifurcation which is a forward bifurcation. Moreover, we showed that model (3) exhibits for some values of the parameters, a varied dynamic, like the permanence and the persistence of the model; a unique positive globally asymptotically stable coexistence equilibrium, connecting with a coexistence of positive saddle equilibrium point; saddle-node equilibria and extinction of two equilibria. We showed by the Bendixson–Dulac criterion that model (3) has no periodic solutions, and then nonexistence of limit cycle is proved. Finally, a numerical simulation is taken to verify the results we obtained.

Li et al. [10] have studied a predator–prey model with nonsmooth switched harvest on the predator, their model presents rich dynamics, whereas the fractional order derivative of the model has an acceleration effect in the dynamics [14]. Naik et al. [12] proved that a discrete-time prey–predator model could exhibit complicated bifurcation phenomena, including period-doubling, Neimark-Sacker and strong resonance bifurcations. Several interesting results (stability of the model, boundedness of the solutions, chaos, etc.) are also obtained by considering three and four dimensional fractional systems [3, 4, 13]. In Refs. [3, 4], the authors proved that the fractional power of the derivative has a significant effect on the dynamic process. Compared with those existing results, our study can be thought as a supplement to model (1). The impact of constant and linear harvesting on system (1) have been extensively studied, but how do nonlinear harvesting affect the dynamic of system (1) is not yet clear. Moreover, nonlinear type harvesting is more realistic and better than the constant and the linear effort harvesting from biological and economic points of view. From the numerical simulations, it may be concluded that low level of harvesting activities leads to the coexistence of the system, but conversely high level of harvesting activities eradicates the entire predator population in short period of time and will lead to the extinction of the system. Illegal fishing and over harvesting on fisheries should be avoided to allow future generations to benefit from marine resources. However, rational initiatives should be conducted to ensure that harvesting levels are at moderate levels in order to protect the entire aquatic ecosystem from collapse. From a practical standpoint, we would like to design harvesting policies to keep the predators and prey from extinction. According to the nature of the system, the decision-makers of a company can develop the best harvesting strategy to ensure the sustainable development of the ecosystem. This research will be useful for further understanding the dynamic complexity of ecosystems when there is the nonlinear type harvesting effect on prey and predator populations.


A part of the bifurcation diagram for model (3) with x2 and rα1 vs α. The solid line presents the curve of the prey with coexistence equilibrium (x2, y2), and the dotted line indicates the curve of the prey with coexistence equilibrium (rα1,0)

Figure 1

A part of the bifurcation diagram for model (3) with x2 and rα1 vs α. The solid line presents the curve of the prey with coexistence equilibrium (x2, y2), and the dotted line indicates the curve of the prey with coexistence equilibrium (rα1,0)

The bifurcation diagram for model (3) with x2, and rα1 vs α, after the forward bifurcation. The solid line presents the curve of the prey with coexistence equilibrium (rα1,0)

Figure 2

The bifurcation diagram for model (3) with x2, and rα1 vs α, after the forward bifurcation. The solid line presents the curve of the prey with coexistence equilibrium (rα1,0)

The phase portrait of system (3) when (x2, y2) is globally asymptotically stable and (rα1,0) is unstable

Figure 3

The phase portrait of system (3) when (x2, y2) is globally asymptotically stable and (rα1,0) is unstable

The phase portrait of system (3) when (x2, y2) meets (rα1,0) at a saddle-node equilibrium

Figure 4

The phase portrait of system (3) when (x2, y2) meets (rα1,0) at a saddle-node equilibrium

The phase portrait of system (3) with a stable equilibrium (rα1,0)

Figure 5

The phase portrait of system (3) with a stable equilibrium (rα1,0)


1.Apima SB, Lawi GO, Kagendo NJ. Density dependent delayed migration for Rosenzweig-Macaurther model with Holling type II predator functional response. Asian Res J Math. 2019; 15(3): 1-12.

2.Das T, Mukherjee RN, Chaudhuri KS. Harvesting of a prey-predator fishery in the presence of toxicity. Appl Math Model. 2009; 33: 2282-92.

3.Daşbaşı B. Stability analysis of an incommensurate fractional-order SIR model. Math Model Numer Simulation Appl. 2021; 1(1): 44-55.

4.Gholami M, Ghaziani RK, Eskandari Z. Three-dimensional fractional system with the stability condition and chaos control. Math Model Numer Simulation Appl. 2022; 2(1): 41-7.

5.Hoekstra J, Van den Bergh JCJM. Harvesting and conservation in a predator-prey system. J Econ Dyn Control. 2005; 29(6): 1097-120.

6.Hu D, Cao H. Stability and bifurcation analysis in a predator-prey system with Michaelis-Menten type predator harvesting. Nonlinear Anal real World Appl. 2017; 33: 58-82.

7.Kar TK. Stability analysis of a prey-predator model incorporating a prey refuge. Commun Nonlinear Sci Numer Simul. 2005; 10: 681-91.

8.Keong AT, Safuan HM, Kavikumar J. The impact of harvesting activities on prey-predator fishery model in the presence of toxin. J Sci Technol. 2018; 10(2): 128-35.

9.Keong AT, Safuan HM, Jacob K. Dynamical behaviors of prey-predator fishery model with harvesting affected by toxic substances. Mathematika. 2018; 34(1): 143-51.

10.Li B, Liu S, Cui J, Li J. A simple predator-prey population model with rich dynamics. Appl Sci. 2016; 6(151): 1-18.

11.Mohdeb N. Stability analysis of prey-predator model with nonlinear harvesting rate. In: Mathematical population: dynamics, ecoepidemiology and evolution MPDEE’21, Centre International de Rencontres Mathmétiques, April 26-29, Luminy; 2021. p. 47. Available from: https://mpdee21.mio.osupytheas.fr/wp-content/uploads/2021/04/Abstract_CMPDE21.pdf

12.Naik PA, Eskandari Z, Yavuz M, Zu J. Complex dynamics of a discrete-time Bazykin–Berezovskaya prey-predator model with a strong Allee effect. J Comput Appl Math. 2022; 413: 1-12: 114401.

13.Naim M, Sabbar Y, Zeb A. Stability characterization of a fractional-order viral system with the non-cytolytic immune assumption. Math Model Numer Simulation Appl. 2022; 2(3): 164-76.

14.Yavuz M, Sene N. Stability analysis and numerical computation of the fractional predator–prey model with the harvesting rate. Fractal and Fractional. 2020; 4(35): 1-22.

15.Berryman AA. The origins and evolution of predator prey theory. Ecology. 1992; 73(5): 1530-5.

16.Brauer F, Soudack AC. Stability regions and transition phenomena for harvested peadator-prey systems. J Math Biol. 1979; 7: 319-37.

17.Beddington JR, Cooke JK. Harvesting from a prey-predator complex. Ecol Model. 1982; 14: 155-77.

18.Berglan H, Wyller J, Burlakov E. Pasture-livestock dynamics with density-dependent harvest and changing environment. Nat resource Model. 2019; 32(4): 1-33: e12213.

19Dumortier F, Libre J, Artés JC. Qualitative theory of planar differential systems. Berlin: Springer; 2006.

20.Hartman P. Ordinary differential equations. New York London Sydney: Wiley; 1964.

21.Meyer KR. The implicit function theorem and analytic differential equations. Lecture Notes Math. 1974; 468: 191-208.

22Layek GC. An introduction to dynamical systems and chaos. New Delhi: Springer; 2015.


The authors like to acknowledge the support of “Direction Générale de la Recherche Scientifique et du Développement Technologique DGRSDT”. MESRS, Algeria.

Corresponding author

Nadia Mohdeb can be contacted at: nadia.mohdeb@univ-bejaia.dz

Related articles