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

Abstract

Purpose

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.

Design/methodology/approach

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.

Findings

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.

Originality/value

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.

Citation

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

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

(1)x.=xr1−xk−axy,y.=caxy−dy.

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 Q_i and F_i are constants and x_i, 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

(2)x.=xr1−xk−axy−αH1x,y.=caxy−dy−βH2y,

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

(3)x.=xr1−xk−axy−αx2,y.=caxy−dy−βy3.

The harvesting function of the prey is quadratic and can be added to the logistic term, and the first equation becomes x.=x(r−xγ)−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

x0>maxrα1,dca.

Then, D=x,y∈R×R; 0<x<x0,0<y<1βcax0−d 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=d−carα1 and Δ=c2a4α12−4β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 d₁ > 0, or in the case d₁ < 0, system (3) admits three positive equilibria: 0,0, rα1,0 and x2,y2 where

x2=rα1−aα1y2 and y2=−ca2α1+Δ2β.

If β=c2a44d1α12, system (3) has two equilibria: 0,0 and rα1,0. When β>c2a44d1α12 and d₁ ≥ 0 there exist two positive equilibria: 0,0 and rα1,0.

Proof. An equilibrium point of system (3) satisfies

(4)xr1−xk−axy−αx2=0,caxy−dy−βy3=0.

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

(5)r1−xk−ay−αx=0

gives x=rα1−aα1y. Solutions of equation

(6)cax−d−βy2=0

are given based on values of d₁ 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 d₁ > 0 or if d₁ < 0; one equilibrium point x1,y1 if d₁ = 0, where

(7)x1=rα1−aα1y1, y1=−ca2α1−Δ2β; 

one equilibrium point x1,−ca22α1β if β=c2a44d1α12; and none if β>c2a44d1α12 and d₁ > 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 d₁ ≥ 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 d₁ < 0, and a locally stable node when d₁ > 0. However, if d₁ = 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 d₁ = 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 W_u tangent at point rα1,0 to the x-axis on which the behavior of system (3) is repulsive.

On one side of W_u, 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 d₁ = 0. First we translate the equilibrium point rα1,0 to the origin by the change of variables X=x−rα1 and Y = y. System (3) becomes

X.=−rX−arα1Y−aXY−α1X2,Y.=caXY−βY3.

By using the transvection u=X+rα1Y, v = Y and by reversing the time τ = −t, we get

(8)u′=ru−a−ca2α1uv+α1u2−ca3α12v2−βaα1v3,v′=−cauv+ca2α1v2+βv3,

where (′) denotes ddt. By exchanging the roles of u and v in system (8) by setting x¯=v, y¯=u, we obtain

(9)x¯′=−cax¯y¯+ca2α1x¯2+βx¯3,y¯′=ry¯−a−ca2α1x¯y¯+α1y¯2−ca3α12x¯2−βaα1x¯3.

We denote Ax,y=−caxy+ca2α1x2+βx3. Equation y¯′=0 gives

(10)α1y¯2+r−a−ca2α1x¯y¯−ca3α12x¯2−βaα1x¯3=0.

The solutions of equation (10) depend on the sign of

Δ-=r2+a+ca2α12x¯2−2ra−ca2α1x¯+4aβx¯3.

Because x¯,y¯ is supposed to be in the neighborhood of the origin, then since x¯ is smallest than r², we obtain Δ->0 and equation (10) has two solutions:

y¯1x¯=12α1−r−a−caα1x¯−Δ-,

and

y¯2x¯=12α1−r−a−caα1x¯+Δ-.

We can show that y¯1 is excluded and y¯=y¯2. By the third-order Taylor development, we get

Ax¯,y¯2x¯=ca2α1x¯2+β−c2a4rα12x¯3+ox¯3.

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 W_u tangent at origin to the y¯-axis such that on W_u, 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 W_u at the point whose ordinate is y¯0, it passes a center manifold W_c tangentially at origin to the x¯-axis. On W_c, 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

x¯′=x¯2+β−c2a4rα12x¯3,y¯′=ry¯,

and C0-conjugate to

x¯′=x¯2,y¯′=y¯.

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 y₂ is positive if and only if d₁ ≤ 0, equilibrium point x2,y2 is not biological meaningful if d₁ > 0; we note that x₂ 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 d₁ ≥ 0, and it is either a linear center, a focus or a node if d₁ < 0; it is locally stable if d₁ < 0 and 1−caα1≥0.

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

Jx,y=−r+ay−arα1−aα1ycaycarα1−aα1y−d−3βy2.

We first assume that d₁ < 0. We have

detJx,y=r−ay−caα1r−ayy+d+3βy2+ca2α1y.

We denote Ix,y=2ca2α1y+d1+3βy2. Its sign depends on that of Δ1≔c2a4α12−3βd1. We have Δ₁ = Δ + βd₁, it is positive. This means that Ix,y has two roots:

y¯=−c2a4α12−Δ23β and y¯¯=−c2a4α12+Δ23β,

where Δ2≔c2a4α12−3βd1. We will compare y₂ with respect to y¯¯. Suppose that y2<y¯¯. Then

(11)2c2a4α12+3−ca2α1+Δ<Δ2.

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βd₁ < 0 which is absurd, and then Ix2,y2>0.

Moreover, we can show that r − ay₂ > 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

trJx,y=r−ay−1+caα1−d−3βy2,

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

If d₁ = 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

h(x,y)=1xy,

for x > 0, y > 0 and denote

f(x,y)=xr1−xk−axy−αx2,g(x,y)=caxy−dy−βy3.

We have

∂(hf)∂xx,y+∂(hg)∂yx,y=−rky−αy−2βyx.

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 d₁ < 0, then the positive equilibrium point x2,y2 is asymptotically stable; it is globally asymptotically stable in D⁺.

2) If d₁ > 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 d₁ < 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 d₁ > 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 d₁ < 0.

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

• A≔x0,y0 a point of the curve y=1βcax−d 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βcax−d;

• 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βcax−d and

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

Let D₁ 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 D₁ is positively invariant compact for system (3). Since in the case d₁ < 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 d₁ ≤ 0 for α ≤ 0.25. The two equilibria (x₂, y₂) 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 d₁ ≥ 0 and β>c2a44d1α12 for 0.25 ≤ α (Figure 2). The positive equilibrium does not exist because y − 2 becomes negative (x₂ is still positive).

If α is less than the value giving d₁ = 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 d₁ = −0.1 < 0. The equilibrium point (x₂, y₂) 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 d₁ = 0. The equilibrium point (x₂, y₂) 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 d₁ = 0.1 > 0. The positive equilibrium point (x₂, y₂) 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.

Figures

Figure 1

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

Figure 2

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

Figure 3

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

Figure 4

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

Figure 5

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