On a predator-prey system interaction under fluctuating water level with nonselective harvesting

Na Zhang; Yonggui Kao; Fengde Chen; Binfeng Xie; Shiyu Li

doi:10.1515/math-2020-0145

Article Open Access

On a predator-prey system interaction under fluctuating water level with nonselective harvesting

Na Zhang , Yonggui Kao , Fengde Chen , Binfeng Xie and Shiyu Li

Published/Copyright: June 4, 2020

Published by

Become an author with De Gruyter Brill

Submit Manuscript Author Information Explore this Subject

$Open Mathematics$

From the journal Open Mathematics Volume 18 Issue 1

Abstract

A predator-prey model interaction under fluctuating water level with non-selective harvesting is proposed and studied in this paper. Sufficient conditions for the permanence of two populations and the extinction of predator population are provided. The non-negative equilibrium points are given, and their stability is studied by using the Jacobian matrix. By constructing a suitable Lyapunov function, sufficient conditions that ensure the global stability of the positive equilibrium are obtained. The bionomic equilibrium and the optimal harvesting policy are also presented. Numerical simulations are carried out to show the feasibility of the main results.

Keywords: predator-prey; permanence; stability; non-selective harvesting; bionomic equilibrium

MSC 2010: 34D20, 34D23, 92D25, 92D40

1 Introduction

Lotka-Volterra systems are one of the most classical and important systems in the field of mathematical biology research and were initially independently proposed in the 1920s by the American biophysicist Lotka when studying a chemical reaction and the Italian mathematician Volterra when studying competition between fish. Generally, we present a kind of predator-prey system between two populations depending on the Lotka-Volterra model as follows:

(1.1) { d x d t = x ( a 10 − a 11 x ) − a 12 x y , d y d t = y ( a 20 − a 22 y ) + a 21 x y ,

where x and y denote the density of the prey population and predator population, respectively, and a i j ( i , j = 1 , 2 … ) are all positive constants. In system (1.1), a 12 x y refers to the quantity of prey that was eaten by the predator per unit time. Thus, a 12 x means the quantity of prey that was eaten by one predator per unit time, which is described as the predation rate and denoted as Φ ( x ) . Then, Φ ( x ) is called the functional response function of prey. Some scholars have devoted their efforts to investigate the predator-prey system in previous studies [1,2,3,4].

Rosenzweig-Macarthur found that the predation rate would be affected not only by the prey population but also by the predator population. In view of this, Rosenzweig-Macarthur improved the Lotka-Volterra model with the consideration of the functional response function as Φ ( x , y ) ,

{ d x d t = f ( x ) − Φ ( x , y ) , d y d t = − e y + k Φ ( x , y ) .

Later, various scholars studied the predator-prey system with different kinds of functional response functions, such as Holling-style functional response function [5,6,7,8,9,10,11,12,13,14,15,16]. Recently, in reference [17], Fellah et al. proposed a predator-prey system with a new kind of functional response function,

(1.2) { d G d t = γ G G ( t ) − min ( r ( t ) G ( t ) B ( t ) + D , γ B ) − m G G 2 ( t ) , d B d t = τ B min ( r ( t ) G ( t ) B ( t ) + D , γ B ) − m B B ( t ) .

This study is based on two interdependent fishes in an artificial lake in French: the pike (B), which is the most important predator, and the roach (G), which is the prey. The functional response function is based on the following general considerations. When a predator attacks a prey, it has access to a certain quantity of food depending on the water level. When the water level is low, during autumn, the predator is more in contact with the prey, and the predation increases. When the water level is high in spring, it is more difficult for the predator to find a prey and the predation decreases. Thus, it is assumed that accessibility function r(t) for the prey is continuous and 1-periodic: the minimum value r ₁ is reached in spring, and the maximum value r ₂ is attained during autumn. The functional response function is given by

min ( r ( t ) G ( t ) B ( t ) + D , γ B ) ,

where D measures other causes of mortality outside of predation and γ B is the maximum predator needs. Obviously, the functional response function given in system (1.2) is reasonable and universal in most of the lakes.

Belkhodja et al. [18] proposed and studied the prey-predator system based on system (1.2) as follows, in which the harvesting of prey population was considered:

(1.3) { d x d t = a x ( t ) ( 1 − x ( t ) K ) min ( b x ( t ) y ( t ) + D , γ ) y ( t ) − q E x ( t ) , d y d t = − d y ( t ) + e min ( b x ( t ) y ( t ) + D , γ ) y ( t ) .

It is different from system (1.2) that the authors provide b = ∫ 0 1 r ( t ) d t as the mean function of predation rate. In reference [18], the authors investigated the dynamics of system (1.3) and established the sufficient criteria for the boundedness, permanence, and predator extinction. Then, the local and the global stability of the equilibrium of system (1.3) are studied. Finally, the authors investigated the bionomic equilibrium and the optimal harvesting policy of system (1.3).

Most of the time, it is difficult to harvest only one kind of fishes independently when we fish in the lakes. As studied in references [19,20,21], we take nonselective harvesting into consideration for the sake of economic profit. We studied a kind of prey-predator system with nonselective harvesting as follows:

(1.4) { d x d t = r x ( 1 − x K ) − min ( b x y + D , γ ) y − q 1 m E x , d y d t = − β y + α min ( b x y + D , γ ) y − q 2 m E y ,

where x(t) and y(t) denote the densities of the prey and predator, respectively, at time t and r , K , b , α , β , γ , q 1 , q 2 , m , D , E are positive constants: r means the intrinsic growth rate of prey, K means the carrying capacity for prey, b is the mean function for the predation rate of prey, β is the death rate of predator, γ means the maximum consumption rate of resource by predator, α means the conversion rate from prey that was eaten by the predator to newborn predators, D measures other causes of mortality outside of predation, q ₁ and q ₂ denote the catchability coefficients of the prey and predator species, respectively, E is the effort devoted to the harvesting of human beings, and m is the fraction of the stock available for harvesting and 0 < m < 1.

This paper is organized as follows. In Section 2 , we obtain the sufficient conditions that ensure the permanence of two populations and the extinction of predator population. In Section 3 , we give all the nonnegative equilibrium points of the system under the necessary assumptions and discuss their local and global stabilities. In Section 4 , we investigated the bionomic equilibrium of the system to guarantee the permanence of two populations under reasonable harvesting efforts. The optimal harvesting policy is discussed in Section 5 . Finally, we conclude with the results obtained in the paper in Section 6 .

2 Permanence

In this section, we study the permanence of two populations. From the standpoint of biology, we are only interested in the dynamics of system (1.4) in ℝ + 2 ≔ { ( x , y ) ∈ ℝ 2 , x ≥ 0 , y ≥ 0 } with the initial condition x ( 0 ) = x 0 > 0 , y ( 0 ) = y 0 > 0 .

Obviously, if E > r q 1 m in system (1.4), then d x d t < 0 , as a result, both prey and predator population will go to extinction. Hence, we assume that

(H0) 0 ≤ E < r q 1 m .

Lemma 1

ℝ + 2 is an invariant set for system (1.4).

For the proof of Lemma 1 refer to the proof of Lemmas 1 and 2 in the literature [18].

Lemma 2

The solution of system (1.4) satisfying the initial value in the ℝ + 2 is bounded and ultimate bound.

Proof

Let w ( t ) = α x ( t ) + y ( t ) . Then the time derivative along the solutions of system (1.4) is given by

d w d t = α d x d t + d y d t = α r x ( 1 − x K ) − α q 1 m E x − β y − q 2 m E y .

Thus,

d w d t + ( β + q 2 m E ) w = α [ − r K x 2 + ( r − q 1 m E + β + q 2 m E ) x ] ≤ α K 4 r ( r − q 1 m E + β + q 2 m E ) 2 ≐ ω .

Let υ = β + q 2 m E . By using the theory of differential inequality, we obtain

0 < w ( t ) ≤ w ( 0 ) exp { − υ t } + ω υ ( 1 − exp { − υ t } ) .

Thus, we get 0 < w ( t ) ≤ max { w ( 0 ) , ω υ } . Moreover, we have lim t → + ∞ w ( t ) ≤ ω υ , from which we can conclude that all the solutions of system (1.4) that start in ℝ + 2 are confined to the region U, U = { ( x , y ) ∈ ℝ + 2 : w ( x , y ) < ω υ + ε } , for any ε.□

For the research of the permanence behavior of the system, we give the definition of permanence and nonpermanence at first.

Definition 1

[18] System (1.4) is called permanent if there exist positive constants 0 < m ⁎ ≤ M ⁎ such that

min { lim inf t → + ∞ x ( t ) , lim inf t → + ∞ y ( t ) } ≥ m ⁎ , max { lim sup t → + ∞ x ( t ) , lim sup t → + ∞ y ( t ) } ≤ M ⁎ ,

for all solutions (x(t), y(t)) of system (1.4) with positive initial values.

The system is said to be nonpermanent if there is a positive solution (x(t), y(t)) of system (1.4) such that

min { lim inf t → + ∞ x ( t ) , lim inf t → + ∞ y ( t ) } = 0 .

Before we establish the persistence for system (1.4), we need to provide a lemma.

Lemma 3

[2] If a, b > 0 and d x d t ≤ ( resp . ≥ ) x ( t ) ( a − b x ( t ) ) with x ₀ > 0, then we have

lim sup t → + ∞ x ( t ) ≤ a b ( resp . lim inf t → + ∞ x ( t ) ≥ a b ) .

The necessary assumptions should be provided before we study the permanence of the populations:

(H1) max { b x 0 y 0 + D , b K ( r − q 1 m E + β + q 2 m E − b ) 2 4 r D ( β + q 2 m E − b ) } < γ ,

(H2) 0 ≤ E < r − b q 1 m ,

(H3) β + q 2 m E r − q 1 m E − b < α b K r D .

Proposition 1

If assumptions (H1)–(H3) hold, then system ( 1.4 ) is permanent, i.e., for any positive solution (x(t), y(t)) with x ( 0 ) > 0 , y ( 0 ) > 0 of system ( 1.4 ), there exist 0 < m i ⁎ ≤ M i ⁎ ( i = 1 , 2 … ) such that

0 < m 1 ⁎ ≤ lim inf t → + ∞ x ( t ) ≤ lim sup t → + ∞ x ( t ) ≤ M 1 ⁎ , 0 < m 2 ⁎ ≤ lim inf t → + ∞ y ( t ) ≤ lim sup t → + ∞ y ( t ) ≤ M 2 ⁎ ,

where

m 1 ⁎ = K r ( r − q 1 m E − b ) , m 2 ⁎ = α b m 1 ⁎ β + q 2 m E − D ,

M 1 ⁎ = K r ( r − q 1 m E ) , M 2 ⁎ = α b M 1 ⁎ β + q 2 m E − D .

Proof

From the first equation of system (1.4), it follows that

d x d t ≤ x [ ( r − q 1 m E ) − r x K ] .

Using Lemma 3, we have

lim sup t → + ∞ x ( t ) ≤ K r ( r − q 1 m E ) ≐ M 1 ⁎ .

Thus, for arbitrary ε 1 > 0 , there exists a positive real number T ₁ such that

(2.1.1) x ( t ) < M 1 ⁎ + ε 1 , for all t ≥ T 1 .

From the second equation of system (1.4) and equation (2.1.1), for all t ≥ T 1 , we have

d y d t ≤ α b ( M 1 ⁎ + ε 1 ) y + D y − ( β + q 2 m E ) y = y y + D [ ( α b ( M 1 ⁎ + ε 1 ) − ( β + q 2 m E ) ( y + D ) ] ≤ y D [ α b ( M 1 ⁎ + ε 1 ) − ( β + q 2 m E ) D − ( β + q 2 m E ) y ] .

Using Lemma 3 and considering the arbitrariness of ε ₁, we obtain

lim sup t → + ∞ y ( t ) ≤ α b M 1 ⁎ β + q 2 m E − D ≐ M 2 ⁎ .

Hence, for arbitrary ε 2 > 0 , there exists a positive real number T 2 > T 1 such that

(2.1.2) y ( t ) < M 2 ⁎ + ε 2 , for all t ≥ T 2 .

Now, we give the following proposition before proving the permanence of system (1.4).□

Proposition 2

If hypothesis (H1) holds, then for all t ≥ 0 we have

b x ( t ) < γ ( y ( t ) + D ) .

Proof

Let

μ ( t ) = b x ( t ) − γ ( y ( t ) + D ) .

From hypothesis (H1), we get μ ( t ) < 0 for all t ≥ 0 . Otherwise, there must exist t 0 > 0 such that μ ( t 0 ) = 0 and d μ d t ( t 0 ) ≥ 0 . From this we get

b x ( t 0 ) y ( t 0 ) + D = γ , or y ( t 0 ) = b x ( t 0 ) γ − D .

Then, we have

d μ d t ( t 0 ) = b d x d t ( t 0 ) − γ d y d t ( t 0 ) = − b r K x 2 ( t 0 ) + b ( r − q 1 m E ) x ( t 0 ) − b γ y ( t 0 ) + β γ y ( t 0 ) − α γ 2 y ( t 0 ) − γ q 2 m E y ( t 0 ) ≤ − b r K x 2 ( t 0 ) + b ( r − q 1 m E ) x ( t 0 ) + γ ( β + q 2 m E − b ) ( b x ( t 0 ) γ − D ) = − b r K [ x ( t 0 ) − K 2 r ( r − q 1 m E + β + q 2 m E − b ) ] 2 + b K 4 r ( r − q 1 m E + β + q 2 m E − b ) 2 − γ D ( β + q 2 m E − b ) .

From hypothesis (H1) we have d μ d t ( t 0 ) < 0 , which leads to a contradiction. Therefore, μ ( t ) < 0 for all t ≥ 0 .□

Therefore, system (1.4) is reduced to the simple form as follows when hypothesis (H1) holds.

(2.1) { d x d t = r x ( 1 − x K ) − b x y y + D − q 1 m E x , d y d t = − β y + α b x y y + D − q 2 m E y .

Now, the proof of the permanence of system (1.4) will continue.

From the first equation of (2.1), we get

d x d t ≥ x ( r − q 1 m E − b − r k x ) .

Using Lemma 3, we get

lim inf t → + ∞ x ( t ) ≥ K r ( r − q 1 m E − b ) ≐ m 1 ⁎ .

Thus, for arbitrary ε 3 > 0 , there exists a positive real number T ₃ such that

(2.1.3) x ( t ) > m 1 ⁎ − ε 3 , for all t ≥ T 3 .

From inequalities (2.1.2) and (2.1.3), we obtain that there exists T 4 = max { T 2 , T 3 } such that

d y d t ≥ y M 2 ⁎ + ε 2 + D [ α b ( m 1 ⁎ − ε 3 ) − ( β + q 2 m E ) D − ( β + q 2 m E ) y ] , for all t ≥ T 4 .

Using Lemma 3 and considering the arbitrariness of ε ₃, we have

lim inf t → + ∞ y ( t ) ≥ α b m 1 ⁎ β + q 2 m E − D ≐ m 2 ⁎ .

Remark 1

From Proposition 1, we obtain, under hypothesis (H1), if

E < min { r − b q 1 m , α b K ( r − b ) − r β D α b K q 1 m + r D q 2 m } ≐ E ⁎ ,

then the populations can be permanent. It means that the harvesting effort of human beings should be smaller than a threshold to ensure that both of the species survive. In example 1, we give a model of system (2.1) with different values for E to study the influence of the harvesting effort on the permanence of the populations. We find that the number of predator populations decreases depending on the increase in the harvesting effort E. The threshold E* we obtained in example 1 is about 7.9955. But from the results (Figure 1) we know, if E ≤ 16.5, the prey and the predator can be permanent (Figure 1(1)–(5)); if E ≥ 16.8, the predator population goes to extinction finally (Figure 1(6)). Obviously, the threshold E* for the permanence of the populations is between 16.5 and 16.8, which is larger than 7.9955. In fact, conditions (H1)–(H3) for the permanence of the system are sufficient but not necessary. Still, the threshold E* is reasonable because the number of predators is less when the harvesting effort E > E*, which is bad for the protection and maintenance of the ecological population.

By denoting

b ⁎ = r D ( β + q 2 m E ) α K ( r − q 1 m E ) ,

we can derive the next proposition easily depending on the result of Proposition 1.

Proposition 3

If b < b*, then the predator population goes to extinction.

Proof

Using the result in Proposition 2, we have

d y d t ≤ y ( − β − q 2 m E + α b M 1 ⁎ y + D ) ≤ y ( − β − q 2 m E + α b M 1 ⁎ D ) .

By the differential inequalities, we get

y ( t ) ≤ y 0 exp { − ( β + q 2 m E ) + α b M 1 ⁎ D } .

The exponent part in the right of the inequalities is negative when b < b*, which means that the predator population y ( t ) → 0 as t → + ∞ . Thus, the predator goes to be extinctive.□

Remark 2

From Proposition 3, we get the following: if b < b*, then the predator population is going to be extinctive. Thus, the predation rate should be larger than a threshold to avoid this case. From a biological point of view, the predator or prey disappears when the predation rate is sufficiently small or large. In example 2, we take different values for b and find that if b ≤ 1.5, then the predator population goes to extinction (Figure 2(1)). But if b ≥ 25.5, from the population ecology perspective, the prey population is almost extinct at some time (Figure 2(6)). If 1.6 ≤ b ≤ 25.5, both of the populations can be permanent (Figure 2(2)–(5)). As we stated in system (1.1), the predation rate depends on the water level in the lake. Thus, if there are too much or too little rain in a period of time, the harvesting of the fishes should decrease or stop for the lake.

$Figure 1 Persistence with different harvesting efforts.$

Figure 1

Persistence with different harvesting efforts.

3 Stability dynamics

Under hypotheses (H0) and (H1), we study the dynamical behavior in this section.

$Figure 2 Persistence with different predation rates.$

Figure 2

Persistence with different predation rates.

3.1 Equilibrium

It is easy to calculate from system (2.1) that there are three equilibria in ℝ + 2 :

The extinction equilibrium point Q 0 = ( 0 , 0 ) .
The trivial equilibrium Q 1 = ( x 1 , 0 ) exists when b < b*, where x 1 = r k ( r − q 1 m E ) .
The positive equilibrium Q ⁎ = ( x ⁎ , y ⁎ ) exists when b > b*, where x* and y* satisfy

(3.1.1) { r − r K x ⁎ − b y ⁎ y ⁎ + D − q 1 m E = 0 , − β + α b x ⁎ y ⁎ + D − q 2 m E = 0 .

From equation (3.1.1), we obtain

x ⁎ = ( β + q 2 m E ) ( y ⁎ + D ) α b ,

and y* satisfies

(3.1.2) A ( y ⁎ ) 2 + B y ⁎ + C = 0 ,

where

A = r ( β + q 2 m E ) > 0 , B = 2 r D ( β + q 2 m E ) − α b K ( r − q 1 m E ) + α b 2 K , C = D [ r D ( β + q 2 m E ) − α b K ( r − q 1 m E ) ] .

If b < b ⁎ , then C > 0 and B > 0 . Thus, there exist two negative real roots or a pair of conjugate imaginary roots with the existence of negative real parts for equation (3.1.2) since A > 0 . Then, only trivial equilibrium Q ¹ exist.

If b > b ⁎ , then C < 0 and system (2.1) has a unique positive equilibrium ( x ⁎ , y ⁎ ) since A > 0 .

After that, we study the dynamical behavior of system (2.1) in terms of the stability of the three equilibria as follows.

3.2 Stability

Theorem 1

(1) The extinction equilibrium point Q 0 ( 0 , 0 ) is a saddle point, then it is unstable.

(2) The trivial equilibrium Q 1 ( x 1 , 0 ) is stable if and only if b < b ⁎ .

(3) The interior equilibrium Q ⁎ ( x ⁎ , y ⁎ ) is locally asymptotically stable when b > b ⁎ .

Proof

The Jacobian matrix of system (2.1) is given by

J ( x , y ) = ( r − 2 r K x − b y y + D − q 1 m E − b D x ( y + D ) 2 α b y y + D − β − q 2 m E + α b D x ( y + D ) 2 ) .

The Jacobian matrix of equilibrium Q 0 ( 0 , 0 ) is given by
(3.2.1) J ( 0 , 0 ) = ( r − q 1 m E 0 0 − β − q 2 m E ) .
The eigenvalues of (3.2.1) are λ 1 0 = r − q 1 m E > 0 and λ 2 0 = − β − q 2 m E < 0 , which means that the equilibrium Q 0 ( 0 , 0 ) is a saddle point and it is unstable.
The Jacobian matrix of equilibrium Q 1 ( x 1 , 0 ) is
(3.2.2) J ( x 1 , 0 ) = ( − ( r − q 1 m E ) − b D x 1 0 − β − q 2 m E + α b K r D ( r − q 1 m E ) ) .
The eigenvalues of (3.2.2) are λ 1 1 = − ( r − q 1 m E ) < 0 and λ 2 1 = − β − q 2 m E + α b K r D ( r − q 1 m E ) . If b < b ⁎ , then λ 2 1 < 0 and the equilibrium Q 1 ( x 1 , 0 ) is stable. On the other hand, if b > b ⁎ , then the equilibrium Q 1 ( x 1 , 0 ) is unstable and the system has a unique positive equilibrium Q ⁎ ( x ⁎ , y ⁎ ) .
The Jacobian matrix of equilibrium Q ⁎ ( x ⁎ , y ⁎ ) is given by

(3.2.3) J ( x ⁎ , y ⁎ ) = ( − r K x ⁎ − b D x ⁎ ( y ⁎ + D ) 2 α b y ⁎ y ⁎ + D − α b x ⁎ y ⁎ ( y ⁎ + D ) 2 ) .

We get

tr J ( x ⁎ , y ⁎ ) < 0 , det J ( x ⁎ , y ⁎ ) > 0 .

Thus, the eigenvalues of (3.2.3), denoted as λ 1 ⁎ and λ 2 ⁎ , are negative or a pair of conjugate complex numbers with negative real parts, which means that the equilibrium Q ⁎ ( x ⁎ , y ⁎ ) is locally asymptotically stable when b > b ⁎ .□

Remark 3

From the discussion of the local stability of the three equilibria we have, the trivial equilibrium Q 0 = ( 0 , 0 ) is unstable and the interior equilibrium Q ⁎ = ( x ⁎ , y ⁎ ) is locally asymptotically stable when b > b ⁎ (Figure 3). But the trivial equilibrium Q 1 = ( x 1 , 0 ) is stable when b < b ⁎ (Figure 4) and unstable when b > b ⁎ . Thus, we believe that there exists a transcritical bifurcation when b = b ⁎ .

Moreover, from the discussion so far, we know that the local stability of the equilibria depends on the relationship between b and b*. In example 5, we simulate the local stability of the system as the value of b changes from 0.5 to 20, where b* = 1.5 (Figure 5). From the result, we can see that the predator goes to extinction and the trivial equilibrium Q ¹ is stable when b ≤ 1.5; both of the prey and the predator can coexist, and the positive equilibrium Q* is locally asymptotically stable when b > 1.5. For the coexistence equilibrium, the number of prey fishes decreases and the number of predator fishes increases first and then decrease when b changes from 2 to 20.

$Figure 3 Stability of the boundary equilibrium.$

Figure 3

Stability of the boundary equilibrium.

$Figure 4 Stability of the coexistence equilibrium.$

Figure 4

Stability of the coexistence equilibrium.

$Figure 5 Stability of system (2.1) depends on b.$

Figure 5

Stability of system (2.1) depends on b.

In fact, we can prove that the positive equilibrium Q ⁎ ( x ⁎ , y ⁎ ) is globally asymptotically stable when it exists by constructing a suitable Lyapunov function.

3.3 Global stability

Theorem 2

The positive equilibrium Q ⁎ ( x ⁎ , y ⁎ ) of system (2.1) is globally asymptotically stable when b > b ⁎ .

Proof

Let

V ( x , y ) = [ ( x − x ⁎ ) − x ⁎ ln ( x x ⁎ ) ] + θ [ ( y − y ⁎ ) − y ⁎ ln ( y y ⁎ ) ] ,

where θ is a positive constant.

Then, V ( x , y ) ≥ 0 and V ( x , y ) = 0 if and only if x = x ⁎ , y = y ⁎ .

The time derivative of V along the solution of system (2.1) is given below:

d V d t = x − x ⁎ x d x d t + θ y − y ⁎ y d y d t = ( x − x ⁎ ) { ( r − q 1 m E ) − r K x − b y y + D − [ ( r − q 1 m E ) − r K x ⁎ − b y ⁎ y ⁎ + D ] } + θ ( y − y ⁎ ) { [ − β − q 2 m E + α b x y + D ] − [ − β − q 2 m E + α b x ⁎ y ⁎ + D ] } = − r k ( x − x ⁎ ) 2 − b D ( x − x ⁎ ) ( y − y ⁎ ) ( y + D ) ( y ⁎ + D ) + α b θ ( y ⁎ + D ) ( x − x ⁎ ) ( y − y ⁎ ) ( y + D ) ( y ⁎ + D ) − α b θ x ⁎ ( y − y ⁎ ) 2 ( y + D ) ( y ⁎ + D ) .

Choosing θ = D α ( y ⁎ + D ) , we obtain

d V d t = − r k ( x − x ⁎ ) 2 − b x ⁎ D ( y + D ) ( y ⁎ + D ) 2 ( y − y ⁎ ) 2 .

Then, we have d V d t ≤ 0 and d V d t = 0 if and only if ( x , y ) = ( x ⁎ , y ⁎ ) . Therefore, the positive equilibrium Q ⁎ ( x ⁎ , y ⁎ ) of system (2.1) is globally asymptotically stable when it exists.□

4 Bionomic equilibrium

The bionomic equilibrium is an amalgamation of the concepts of biological equilibrium and economic equilibrium. The economic equilibrium is said to be obtained when the TR (the total revenue obtained by selling the harvested biomass) equals Tc (the total cost for the effort devoted to harvesting).

Let c = constant fishing cost per unit effort, p ₁ = the price per unit biomass of prey fish, p ₂ = the price per unit biomass of predator fish.

Then, the net revenue at any time is given by

π ( x , y , E ) = p 1 q 1 m E x + p 2 q 2 m E y − c E ,

where E ∈ [ 0 , E max ] and E max means the maximum harvesting effort that humans provide.

Then, the bionomic equilibrium Q ∞ ( x ∞ , y ∞ , E ∞ ) is given as follows:

(4.1.1) r ( 1 − x K ) − b y y + D − q 1 m E = 0 − β + α b x y + D − q 2 m E = 0 p 1 q 1 m E x + p 2 q 2 m E y − c E = 0 }

We obtain that

x ∞ = c − p 2 q 2 m 2 y ∞ p 1 q 1 m 1 , E ∞ = 1 q 2 m 2 ( α b x ∞ y ∞ + D − β ) ,

where y _∞ satisfy

(4.1.2) A 1 y 2 + B 1 y + C 1 = 0 ,

where

A 1 = r p 2 q 2 2 m > 0 , B 1 = [ q 1 m K ( β p 1 q 1 + r p 1 q 2 ) − q 2 r c ] − p 1 q 1 q 2 m b K + α p 2 q 2 m b K q 1 + r D p 2 q 2 2 m , C 1 = D [ q 1 m K ( β p 1 q 1 + r p 1 q 2 ) − q 2 r c ] − α b c K q 1 .

It is easy to check that if C 1 < 0 , then B 1 2 − 4 A 1 C 1 > 0 since A 1 > 0 , there can be a unique positive solution y _∞. Moreover, if

(H4) y ∞ < c p 2 q 2 m ,

there is a unique bionomic equilibrium of (4.1.1). The sufficient and necessary condition for C 1 < 0 is given by

(H5) p 1 q 1 m D K ( β q 1 + r q 2 ) α b K q 1 + r D q 2 < c .

Theorem 3

If b > b ⁎ , and conditions (H4) and (H5) hold, then the bionomic equilibrium exists.

Remark 4

For the predator-prey system with nonselective harvesting we studied, the existence of the bionomic equilibrium is of great significance. If E = E ∞ , then the net income from fishing is 0. Thus, the fishermen have to cut back on their fishing efforts to make profits, which means E < E ∞ . For an open access fishery, the situation will not last long because more and more fishermen are attracted by interests and join in the fishing. Hence, the harvesting effort E increases and exceeds E ∞ for some time, and then the total cost utilized in harvesting the prey fish and predator fish would outstrip the total revenues obtained from the fishing. Obviously, some of the fishermen would give up fishing because they cannot afford the loss. Therefore, the fishing effort E decreases and E > E ∞ cannot be maintained indefinitely. Briefly, when E < E ∞ or E > E ∞ , the fishing effort will change in the direction of E = E ∞ . Therefore, the existence of economic equilibrium provides a yardstick for the harvesting effort. To verify the existence of the bionomic equilibrium, example 6 (Figure 6) is given in appendix. We use computer tools to obtain a high precision solution rather than an exact solution.

Furthermore, as we can see in the market economy, the changes in fishing costs and fish prices have an impact on the fishery. For this reason, we conducted the following research and obtained three tables.

$Figure 6 Bionomic equilibrium.$

Figure 6

Bionomic equilibrium.

Take r = 12 , K = 30 , b = 8 , α = 0.9 , β = 2 , q 1 = q 2 = 0.6 , m = 0.8 , D = 5 .

From Tables 1–3 we can see that the harvesting effort E _∞ in a bionomic equilibrium drops off with the increase in the fishing cost. But E _∞ increases gradually when the price of prey fish or predator fish increases. Thus, commercial fishing will be reduced with the increase in the fishing cost or the low price of the fishes. On the contrary, E _∞ will decrease.

Table 1

The fishing cost c changes from 28 to 70 with p 1 = 7 , p 2 = 12

c	28	30	40	50	60	70
x _∞	7.7216	7.8122	8.2376	8.6214	8.9684	9.2833
y _∞	0.3568	0.6512	2.1392	3.6514	5.1851	6.7375
E _∞	17.4551	16.5693	13.1414	10.7812	9.0415	7.6970

Table 2

The price of prey fish p ₁ changes from 5 to 11 with c = 60 , p 2 = 12

p ₁	5	6	7	8	9	10	11
x _∞	9.2495	9.1098	8.9684	8.8256	8.6816	8.5368	8.3915
y _∞	6.5627	5.8617	5.1851	4.5329	3.9055	3.3027	2.7245
E _∞	7.8325	8.4140	9.0415	9.7203	10.4562	11.2563	12.1287

Table 3

The price of predator fish p ₂ changes from 10 to 16 with c = 60 , p 1 = 7

p ₂	10	11	12	13	14	15	16
x _∞	9.1562	9.0561	8.9684	8.8910	8.8221	8.7604	8.7047
y _∞	6.0907	5.6007	5.1851	4.8279	4.5175	4.2452	4.0042
E _∞	8.2170	8.6477	9.0415	9.4034	9.7374	10.0468	10.3344

5 Optimal harvesting policy

In order to obtain the optimal harvesting policy [22] of system (2.1), we consider the present value J of a continuous time stream of revenues

J = ∫ 0 + ∞ exp { − δ t } ( p 1 q 1 m E x + p 2 q 2 m E y − c E ) d t ,

where δ is the instantaneous annual rate of discount and E ∈ [ 0 , E max ] , where E _max means the maximum harvesting effort that humans provide.

In this part, we study the optimal harvesting control of system (2.1), which is denoted as E _δ and the corresponding states of the populations are x _δ and y _δ.

We construct the Hamiltonian function as follows:

H = exp { − δ t } ( p 1 q 1 m E x + p 2 q 2 E m y − c E ) + λ 1 [ r x ( 1 − x K ) − b x y y + D − q 1 m E x ] + λ 2 ( − β y + α b x y y + D − q 2 m E y ) ,

where λ ₁ and λ ₂ are the adjoint variables.

If E = 0, it is not meaningful to discuss the optimal harvesting policy. Today, with the highly developed technology, we can make sufficient harvesting efforts so that the population go to extinction. Hence, by the maximum principle, the maximum must occur at ( x δ , y δ , E δ ) , where E δ ∈ ( 0 , E max ) and they satisfy

(5.1) { ∂ H ∂ E = exp { − δ t } ( p 1 q 1 m x + p 2 q 2 m y − c ) − λ 1 q 1 m x − λ 2 q 2 m y = 0 , ∂ H ∂ x = exp { − δ t } p 1 q 1 m E + λ 1 ( r − 2 r K x − b y y + D − q 1 m E ) + λ 2 α b y y + D = 0 , ∂ H ∂ y = exp { − δ t } p 2 q 2 m E − λ 1 b D x ( y + D ) 2 + λ 2 ( − β + α b D x ( y + D ) 2 − q 2 m E ) = 0 , ∂ H ∂ λ 1 = r x ( 1 − x K ) − b x y y + D − q 1 m E x = 0 , ∂ H ∂ λ 2 = − β y + α b x y y + D − q 2 m E y = 0 .

For such a nonlinear homogeneous equation of higher order as above, it is difficult to provide an expression of its solution. Taking into consideration another view point, we can obtain the optimal harvesting control by using least squares method with computer tools, such as MATLAB. Take r = 12 , K = 30 , b = 8 , α = 0.9 , β = 2 , q 1 = q 2 = 0.6 , m = 0.8 , D = 5 , p 1 = 7 , p 2 = 12 , c = 60 , for example. We obtain the optimal harvesting control as E δ = 3.1564 and ( x δ , y δ ) = ( 10.7536 , 17.0249 ) (Figure 7).

$Figure 7 Optimal harvesting policy.$

Figure 7

Optimal harvesting policy.

6 Conclusions

A brief discussion is presented to conclude this work. We study a prey-predator system interaction under fluctuating water level with nonselective harvesting. The permanence of two populations, the local stability, and the global stability of the equilibria are investigated for the study of the behavioral dynamics of the system. The bionomic equilibrium and the optimal harvesting policy are discussed. After each conclusion, the corresponding simulation was given to prove that our conclusion is correct.

Through our research, we find that the predation rate plays a decisive role in both population persistence and system stability. The predation rate depends on the water level in the lake. The predation rate decreases when the water level increases and increases when the water level decreases. Thus, the fishermen should adjust their catch of fish according to seasonal changes for the sake of long-term survival of fish stocks.

In addition, the harvesting effort of human beings has a more significant impact on population persistence. Overfishing can lead to rapid population reduction or even extinction. For example, in South China Sea, there has been overfishing of fishery resources in recent years. As a result, the fishery resources in the South China Sea have declined seriously, and the main economic fish resources have been greatly reduced. There are serious problems in fisheries, such as the decline of fishing efficiency, the closure of fishing vessels, and the decline of fishermen’s income. In order to protect fishery resources and the long-term interests of fishermen in the South China Sea, the government and the South China Sea fishery administration have taken various measures to protect fishery resources, such as establishing prohibited fishing areas, suspending fishing periods, and so on.

7 Numerical simulations

Here, some numerical examples are presented to illustrate the practicability of the theoretical analysis provided in our paper.

Example 1

Take r = 12 , K = 30 , b = 6 , D = 5 , α = 0.9 , β = 2 , q 1 = q 2 = 0.6 , m = 0.8 and the initial condition is ( x ( 0 ) , y ( 0 ) ) = ( 25 , 15 ) (Figure 1).

Example 2

Take r = 12 , K = 25 , D = 5 , E = 4 , α = 0.4 , β = 2 , q 1 = q 2 = 0.25 , m = 0.8 and the initial condition is ( x ( 0 ) , y ( 0 ) ) = ( 15 , 25 ) (Figure 2).

Example 3

Take r = 12 , K = 15 , b = 4 , D = 5 , E = 5 , α = 0.25 , β = 2 , q 1 = q 2 = 0.2 , m = 0.8 . The initial conditions are (x(0), y(0)) = (30, 25) (30, 20), (28, 15), (25, 13), (20, 15), (10, 25), (6, 10), (8, 20), (6, 6), and (7, 15) (Figure 3).

Example 4

Take r = 12 , K = 30 , b = 6 , D = 5 , E = 5 , α = 0.9 , β = 2 , q 1 = q 2 = 0.6 , m = 0.8 . The initial conditions are (x(0), y(0)) = (30, 25) (30, 20), (28, 15), (25, 13), (20, 15), (10, 25), (6, 10), (8, 20), (6, 6), and (7, 15) (Figure 4).

Example 5

Take r = 12 , K = 25 , D = 5 , E = 4 , α = 0.4 , β = 2 , q 1 = q 2 = 0.25 , m = 0.8 , and the initial condition is ( x ( 0 ) , y ( 0 ) ) = ( 30 , 25 ) . The values of b for the lines from right to left are 0.5, 1, 1.5, 2, 4, 8, 12, and 20 respectively (Figure 5).

Example 6

Take r = 12 , K = 30 , b = 8 , D = 5 , α = 0.6 , β = 2 , q 1 = 0.6 , q 2 = 0.5 , m = 0.8 , p 1 = 7 , p 2 = 12 , c = 60 (Figure 6).

Acknowledgments

The authors gratefully acknowledge the anonymous referee’s comments and patient review. This work was supported by the Jiangxi Provincial Education Department funding GJJ170567, Jiangxi Provincial Education Department funding GJJ170566, National Science Foundation of China (No. 11561028 and 11801238) and the Innovation and Entrepreneurship Training Program of Jiangxi University of Technology (DC2018-071).

References

[1] Sze Bi Hsu and Tzy-Wei Huang, Global stability for a class of predator-prey systems, SIAM J. Appl. Math. 55 (1995), no. 3, 763–783, 10.1016/j.cam.2004.10.001.Search in Google Scholar

[2] Fengde Chen, On a nonlinear nonautonomous predator-prey model with diffusion and distributed delay, J. Comput. Appl. Math. 180 (2005), no. 1, 33–49, 10.1016/j.cam.2004.10.001.Search in Google Scholar

[3] Liujuan Chen, Fengde Chen, and Yiqin Wang, Influence of predator mutual interference and prey refuge on Lotka-Volterra predator-prey dynamics, Commun. Nonlinear Sci. Numer. Simul. 18 (2013), no. 11, 3174–3180, 10.1016/j.cnsns.2013.04.004.Search in Google Scholar

[4] A. Moussaoui and S. M. Bouguima, A prey-predator interaction under fluctuating level water, Math. Methods Appl. Sci. 38 (2015), no. 1, 123–137, 10.1002/mma.3055.Search in Google Scholar

[5] Jinhuang Chen and Runxin Wu, A commensal symbiosis model with non-monotonic functional response, Commun. Math. Biol. Neurosci. 2017 (2017), 5.Search in Google Scholar

[6] John B. Collings, The effects of the functional response on the bifurcation behavior of a mite predator-prey interaction model, J. Math. Biol. 36 (1997), no. 2, 149–168, 10.1007/s002850050095.Search in Google Scholar

[7] Xianning Liu and Lansun Chen, Complex dynamics of holling type ii Lotka-Volterra predator-prey system with impulsive perturbations on the predator, Chaos Solitons Fractals 16 (2003), no. 2, 311–320.10.1016/S0960-0779(02)00408-3Search in Google Scholar

[8] T. K. Kar and H. Matsuda, Global dynamics and controllability of a harvested prey-predator system with holling type iii functional response, Nonlinear Anal. Hybrid Syst. 1 (2007), no. 1, 59–67, 10.1016/j.nahs.2006.03.002.Search in Google Scholar

[9] Fengying Wei, Existence of multiple positive periodic solutions to a periodic predator-prey system with harvesting terms and holling iii type functional response, Commun. Nonlinear Sci. Numer. Simul. 16 (2011), no. 4, 2130–2138, 10.1016/j.cnsns.2010.08.028.Search in Google Scholar

[10] J. H. P. Dawes and M. O. Souza, A derivation of holling’s type i, ii and iii functional responses in predator-prey systems, J. Theoret. Biol. 327 (2013), 11–22, 10.1016/j.jtbi.2013.02.017.Search in Google Scholar PubMed

[11] Jun Zhou and Chan-Gyun Kim, Positive solutions for a Lotka-Volterra prey-predator model with cross-diffusion and holling type-ii functional response, Sci. China Math. 57 (2014), no. 5, 991–1010, 10.1007/s11425.013.4711.0.Search in Google Scholar

[12] Zhiqi Lu and Xia Liu, Analysis of a predator-prey model with modified holling-tanner functional response and time delay, Nonlinear Anal. Real World Appl. 9 (2008), no. 2, 641–650, 10.1016/j.nonrwa.2006.12.016.Search in Google Scholar

[13] Jiangbin Chen, Bifurcations of degenerate singular point of a predator-prey system with holling-iv functional response, Ann. Differ. Equations 20 (2014), no. 4, 386–393, 10.1016/s0960.0779(02)00408.3.Search in Google Scholar

[14] Yansen Lv and Zengji Du, Existence and global attractivity of a positive periodic solution to a Lotka-Volterra model with mutual interference and holling iii type functional response, Nonlinear Anal. Real World Appl. 12 (2011), no. 6, 3654–3664, 10.1016/j.nonrwa.2011.06.022.Search in Google Scholar

[15] Hafizul Molla, Md Sabiar Rahman, and Sahabuddin Sarwardi, Dynamics of a predator-prey model with holling type ii functional response incorporating a prey refuge depending on both the species, Int. J. Nonlinear Sci. Numer. Simul. 20 (2019), no. 1, 1–16, 10.1515/ijnsns.2017.224.Search in Google Scholar

[16] Runxin Wu, Lin Li, and Xiaoyan Zhou, A commensal symbiosis model with holling type functional response, J. Math. Comput. Sci. 16 (2016), no. 3, 364–371, 10.22436/jmcs.016.03.06.Search in Google Scholar

[17] N. Chiboub Fellah, S. M. Bouguima, and A. Moussaoui, The effect of water level in a prey-predator interactions: a nonlinear analysis study, Chaos Solitons Fractals 45 (2012), no. 3, 205–212, 10.1016/j.chaos.2011.12.002.Search in Google Scholar

[18] K. Belkhodja, A. Moussaoui, and M. A. Aziz Alaoui, Optimal harvesting and stability for a prey-predator model, Nonlinear Anal. Real World Appl. 39 (2018), 321–336, 10.1016/j.nonrwa.2017.07.004.Search in Google Scholar

[19] Chaoquan Lei, Dynamic behaviors of a non-selective harvesting may cooperative system incorporating partial closure for the populations, Commun. Math. Biol. Neurosci. 2018 (2018), 12.10.31193/ssap.isbn.9787520118873Search in Google Scholar

[20] Kunal Chakraborty, Sanjoy Das, and T. K. Kar, On non-selective harvesting of a multispecies fishery incorporating partial closure for the populations, Appl. Math. Comput. 221 (2013), 581–597, 10.1016/j.amc.2013.06.065.Search in Google Scholar

[21] Qifa Lin, Dynamic behaviors of a commensal symbiosis model with non-monotonic functional response and non-selective harvesting in a partial closure, Commun. Math. Biol. Neurosci. 2018 (2018), 4.Search in Google Scholar

[22] Colin Whitcomb Clark, Mathematical Bioeconomics: The Optimal Management of Renewable Resources (2nd ed.), John Wiley and Sons, New York, 1990.Search in Google Scholar

Received: 2019-03-30

Revised: 2020-02-29

Accepted: 2020-03-20

Published Online: 2020-06-04

This work is licensed under the Creative Commons Attribution 4.0 International License.

Articles in the same Issue

https://doi.org/10.1515/math-2020-0145

Keywords for this article

predator-prey; permanence; stability; non-selective harvesting; bionomic equilibrium

Creative Commons

BY 4.0