key: cord-0883547-fooljxoe
authors: Kaviya, R.; Muthukumar, P.
title: Dynamical analysis and optimal harvesting of conformable fractional prey–predator system with predator immigration
date: 2021-05-16
journal: Eur Phys J Plus
DOI: 10.1140/epjp/s13360-021-01559-w
sha: 99c79a2f30e0a8118099df048745fa67396ca30e
doc_id: 883547
cord_uid: fooljxoe

Aim of this work is to study the four species various fractional-order prey–predator or Lotka–Volterra (LV) system with both immigration and harvesting effects. The existence and uniqueness, uniform boundedness, persistence, permanence, and extinction of this system solution are analyzed. The stability behavior of the system is obtained with the help of the Routh–Hurwitz (RH) stability criterion. The small changes in fractional-order values can produce a significant impact on the stability of the system is confirmed. This work verifies that the small amount of immigration effect can change the dynamic nature of the LV system. Numerical results are given to illustrate the obtained theoretical results of the stability analysis. The bionomic equilibrium points of the system are attained with their feasibility conditions. To get the optimal amount of harvesting effect with the Pontryagin’s maximum principle, the harvesting parameter is considered as the control parameter.

The ecosystem has a wide diversity in the population-community. Such a populationcommunity in an ecosystem is unlikely to survive without dependence on other communities. In that case, there must be a significant relationship between different populationcommunities. Interactions between different populations are generally characterized by mutualism, competition, parasitism, and predation. In the above types of relationships, predatory activity is widespread in real-time. The predatory relationship between predator and prey is widely seen in real-life such as Lotka and Volterra (LV) system (see [1] ). This LV system is employed to characterize the dynamical behavior of the species by their predatory relationship (see [2, 3] ). Recently, mathematicians and ecologists have studied a wide variety of prey-predator systems (see [4] [5] [6] [7] [8] ).

Many researchers studied a large number of fractional-order aspects of the biological process (see [9] [10] [11] [12] [13] ). Those aspects should not be shown and explained by the integer-order calculus, but the fractional calculus can do it. The fractional-order derivative is related to a e-mail: r.kaviya26@gmail.com (corresponding author) b e-mail: pmuthukumargri@gmail.com both immigration and harvesting effects. With this motive, this article investigates the new problem of a conformable fractional-order LV system with the immigration and harvesting effects. The main contributions of this work are given as follows:

• The proposed conformable fractional-order four species LV system includes prey, secondary-level predator, third-level predator, and top-level predator with the various fractional-orders (α 1 , α 2 , α 3 , α 4 ) ∈ (0, 1]. • In this model, each species having intercommunication with another three species. Moreover, the immigration effect is applied to the predator species and harvesting effects are implemented in all species. • The impact of the fractional orders and immigration effect on the dynamical behavior of the LV system is investigated. From this investigation, we can see how the small changes in fractional orders and immigration effect can change the dynamical behavior of the system. • In this study, we have newly achieved the optimal harvesting of the fractional-order LV system by the Pontryagin's maximum principle at the bionomic equilibrium points.

This paper is organized as follows: Sect. 2 gives the preliminary results based on the conformable fractional-order derivative and constructs the four species conformable fractionalorder LV system with immigration and harvesting effects. Section 3 is dedicated to analyzing the existence and uniqueness of the solution of the system and described the solution behavior, uniform boundedness, persistence, permanence, and extinction results. Section 4 verifying the impact of fractional orders and the immigration on stability behavior of the system with the numerical example. In Sect. 5, the optimal harvesting policy of the LV system is obtained by the Pontryagin's maximum principle at the non-trivial bionomic equilibrium points.

In this section, the required basic definition and theorem are defined and a description of the proposed conformable fractional-order LV system is presented.

Definition 1 (See [15] ) Given a differentiable function f :

Theorem 1 (See [15] ) Let α ∈ (0, 1] and f be differentiable at a point t

Let x 1 (t), x 2 (t), x 3 (t) and x 4 (t) denote the densities of prey, secondary-level predator, thirdlevel predator, and top-level predator, respectively, at time t ∈ (t 0 , ∞), t 0 ≥ 0. For the notation purpose, x j (t) is simply notated as x j , j = 1, 2, 3, 4. Assume that x 1 is a prey of its predators x 2 , x 3 and x 4 ; x 2 is a prey of its predators x 3 and x 4 ; and x 3 is a prey of its predator x 4 . Species x j s grow logistically with the carrying capacity k j , j = 1, 2, 3, 4, respectively. The immigration effect I j is applied to the predator species x j , j = 2, 3, 4, respectively. The harvesting effect q j e j is applied to the species x j , j = 1, 2, 3, 4, respectively. With the above assumptions, the following four-dimensional conformable fractional-order LV system is configured,

Logistic growth

;

Logistic growth 

with the initial condition x j (t 0 ) > 0 at t 0 ≥ 0, for all j = 1, 2, 3, 4. The non-negative parameters r j > 0 and k j > 0 are intrinsic growth rate and carrying capacity of each x j , j = 1, 2, 3, 4, respectively, for simplification r j k j is denoted as a j j and a ji > 0 denotes the interaction coefficient between the each species x j , i = j, i, j = 1, 2, 3, 4. The parameter I j > 0 denotes the immigration rate of the each predator species x j , j = 2, 3, 4. Also, the parameters d j > 0, q j > 0 and e j > 0 denote the death rate, catch-ability coefficient, and total effect applied for harvesting of the each species x j , j = 1, 2, 3, 4, respectively.

For all j = 1, 2, 3, 4, the LV system (1) of equations can be defined in the following general form:

Here, f j :

e j , i, j = 1, 2, 3, 4, and I 1 = 0. According to Theorem 1 in [15] , for t ∈ (t 0 , ∞), t 0 ≥ 0, j = 1, 2, 3, 4, the fractionalorder LV system (2) can be written as,

Here, l j = t α j −1 , j = 1, 2, 3, 4. The integer-order LV system (3) (corresponding to the fractional-order LV system (1)) can be written in the general form with the initial condition

Here,

Remark 1 Comparing the LV systems (2) and (4), we can understand that the integer-order LV system (4) is the equivalent form of the fractional-order LV system (2).

This section verifying the existence and uniqueness of the system solution, uniform boundedness, persistence, permanence, and extinction of solution of the LV system (2). (4) is Lipschitz continuous and there exists a constant k ∈ (0, 1) sufficient for the existence and uniqueness of solution

The solution X (t) of the LV system (2) with the initial condition X (t 0 ) and can be written as,

In Eq.

Here, T 0 ∈ (t 0 , t), and T denote the transpose of the matrix. For all i, j = 1, 2, 3, 4,

Choosing supremum norm, for continuous function

Then, we can easily obtained that,

Here, k = T 0 max{ζ 1 , ζ 2 , ζ 3 , ζ 4 } and ζ j = l j (r j + I j − d j − q j e j + A(2a j j + i = j a i j + i = j a ji )), i, j = 1, 2, 3, 4. From the inequality (6), the mapping η is a contraction if k ∈ (0, 1). By contraction principle, η has a unique fixed point X (t) which satisfies (5) . Hence, X (t) is unique solution of the LV system (2).

Theorem 3 (Uniformly bounded) If the conditions r j + I j > d j + q j e j , j = 1, 2, 3, 4, and a i j = a ji , for i > j, i, j = 1, 2, 3, 4 hold, then the LV system (2) is uniformly bounded.

Here, μ j = r j + I j − d j − q j e j , j = 1, 2, 3, 4. By applying the Gronwall's inequality, then

Hence, all the species are uniformly bounded in R 4 + for any initial X (t 0 ) > 0.

Remark 2 Biologically, it is confirmed that the LV system (2) is uniformly bounded when the following two conditions are satisfied.

(i) Sum of the growth rate and immigration rate is equal to the sum of death rate and harvesting effect. (ii) The prey consumption rates are equal to the corresponding predation rates.

The LV system (2) is said to be weakly persistent if every solution X (t) satisfies, the following two conditions are satisfied.

The LV system (2) is said to be strongly persistent if every solution X (t) satisfies the following two conditions.

Theorem 4 (Persistent and permanent) If M j > m j > 0, j = 1, 2, 3, 4 and 0 < x i ≤ 1, i = 2, 3, 4 hold, then the LV system (2) is persistent and permanent, i.e., any positive solution X (t) of the LV system (2) 

Proof As the solution X (t) of the LV system (2) is positive, it follows that,

By the Lemma 3 in [27] , and for arbitrary ε j > 0, there exist positive constants

According to the assumptions, m j > 0, j = 1, 2, 3, 4, and 0 < x i ≤ 1, i = 2, 3, 4, the LV system (2) is equivalently written as,

Here, for all i, j = 1, 2, 3, 4,

Then, M j > m j , j = 1, 2, 3, 4. By the Lemma 3 in [27] , and for arbitrary δ j > 0, there exist positive constants

From (7) and (8), it is clear that the LV system (2) is persistent and permanent.

Remark 3 From this analysis, we have observed that, if the predator species' (x 2 , x 3 , x 4 ) level should not exceed 1, then the LV system (2) is both persistent and permanent.

The LV system (2) is said to be extinction if there is a positive solution X (t) that satisfies min j {lim t→∞ x j (t)} = 0, j = 1, 2, 3, 4.

If the following two conditions, (i) x i ≤ 1 (ii) r j + I j + i< j a ji < d j + q j e j hold for all j = 1, 2, 3, 4, i = 1, 2, 3, then the species x j will be extinct as t → ∞.

Proof As the solution X (t) of the LV system (2) is positive, it follows that:

Since, r 1 < d 1 + q 1 e 1 , we get lim t→∞ x 1 (t) = 0. It means that the species x 1 is extinct as t → ∞. Using the assumption, x i ≤ 1, i = 1, 2, 3, then the LV system (2) can be written as,

Using the assumption, r j + I j + i< j a ji < d j + q j e j , for j = 1, 2, 3, 4, i = 1, 2, 3, then lim t→∞ x j (t) = 0, j = 2, 3, 4. Hence, the species x j , j = 2, 3, 4 are extinct as t → ∞. Remark 4 As t → ∞ the prey species x 1 will be extinct if sum of the death rate and harvesting rate of x 1 is greater the growth rate of x 1 . As t → ∞, the other three predator species x j , j = 2, 3, 4 are extinct if their corresponding prey levels x i , i < j should not exceed 1 and the sum of the death rate and harvesting rate of x j is greater than the sum of the growth, immigration and consumption rates of x j , j = 2, 3, 4.

This section analyzes the stability properties of the LV system (2) at the positive interior equilibrium point

Here, the stability analysis of the LV system (2) is obtained via the Routh-Hurwitz stability test.

First to construct the Jacobian matrix G = (g ji ) for all i = j and g j j = 0 of the LV system (2) at the positive interior equilibrium point

for analyzing the stability. G = ⎡ ⎢ ⎢ ⎣ g 11 g 12 g 13 g 14 g 21 g 22 g 23 g 24 g 31 g 32 g 33 g 34 g 41 g 42 g 43 g 44

For, i, j = 1, 2, 3, 4,

The characteristic equation of the Jacobian matrix G at the positive interior equilibrium point

Here, G 1 = −(g 11 + g 22 + g 44 ); G 2 = g 11 (g 22 + g 44 ) + g 44 (g 22 + g 33 ) + x * 1 l 1 (a 13 a 31 x * 3 l 3 + a 14 a 41 x * 4 l 4 ) + x * 2 l 2 (a 24 a 42 x * 4 l 4 + a 12 a 21 x * 1 l 1 ) + x * 3 l 3 (a 34 a 43 x * 4 l 4 + a 23 a 32 x * 2 l 2 ); G 3 = −[g 44 (g 11 ( g 22 + g 33 ) + g 22 g 33 ) + g 11 (a 34 a 43 x * 3 x * 4 l 3 l 4 ) + g 11 (a 23 a 32 x * 3 l 2 l 3 + a 24 a 42 x * 2 x * 4 l 2 l 4 )x * 2 + g 22 (a 34 a 43 x * 3 x * 4 l 3 l 4 ) + a 23 a 32 x * 2 x * 3 l 2 l 3 g 44 + a 23 a 34 a 42 x * 2 x * 3 x * 4 l 2 l 3 l 4 − a 24 a 32 a 43 x * 2 x * 3 x * 4 l 2 l 3 l 4 + a 24 a 42 x * 2 x * 4 l 2 l 4 g 33 + a 12 a 21 x * 1 x * 2 l 1 l 2 g 44 + a 12 a 23 a 31 x * 1 x * 2 x * 3 l 1 l 2 l 3 + a 12 a 24 a 41 x * 1 x * 2 x * 4 l 1 l 2 l 4 − a 13 a 21 a 32 x * 1 x * 2 x * 3 l 1 l 2 l 3 + a 13 a 31 x * 1 x * 3 l 1 l 3 g 44 + a 13 a 34 a 41 x * 1 x * 3 x * 4 l 1 l 3 l 4 − a 14 a 21 a 42 x * 1 x * 2 x * 4 l 1 l 2 l 4 − a 14 a 31 a 43 x * 1 x * 3 x * 4 l 1 l 3 l 4 + a 14 a 41 x * 1 x * 4 l 1 l 4 g 33 ]; G 4 = g 11 g 22 g 33 g 44 + g 11 g 22 (a 34 a 43 x * 3 x * 4 l 3 l 4 ) + g 11 (a 23 a 32 x * 2 x * 3 l 2 l 3 g 44 + a 23 a 34 a 42 x * 2 x * 3 x * 4 l 2 l 3 l 4 − a 24 a 32 a 43 x * 2 x * 3 x * 4 l 2 l 3 l 4 + a 24 a 42 x * 2 x * 4 l 2 l 4 g 33 ) + a 12 a 21 x * 1 x * 2 l 1 l 2 g 33 g 44 + a 12 a 23 a 31 x * 1 x * 2 x * 3 l 1 l 2 l 3 g 44 12 a 23 a 34 a 41 − a 12 a 24 a 31 a 43  + a 12 a 21 a 34 a 43 − a 13 a 21 a 34 a 42 + a 13 a 24   × (a 31 a 42 − a 34 a 41 ) + a 14 a 21 a 32 a 43 − a 14 a 23 (a 31 a 42 − a 32 a 41 ) )

Thus, the above discussion gives the following stability theorem.

Remark 5 From Theorem 6, it can be observed that the fractional-orders α j (here, α j viewed in terms of l j = t α j −1 , j = 1, 2, 3, 4) also perform an essential role in stability analysis. It is proved that the stability of the LV system (2) depends on fractional-orders of the LV system (2). In Sect. 4.2, stability behavior of the LV system (2) with the various fractional orders is graphically illustrated.

If the condition r j + I j + i< j (a ji x i ) ≥ +a j j x j + i> j (a ji x i )+d j +q j e j , j = 1, 2, 3, 4, hold for the LV system (2), then the positive equilibrium point

To demonstrate the impact of fractional-orders in stability of the LV system (2) with the parameter set 1 and a different set of fractional orders. Parameter set-1 Consider the following non-negative parameters (r 1 , r 2 , r 3 , r 4 ) = (0.05, 0.009, 0.005, 0.05),  (k 1 , k 2 , k 3 , k 4 ) = (0.75, 0.7, 0.5, 0.2), (a 12 , a 13 , a 14 ) = (0.01, 0.5, 0.01), (a 21 , a 23 , a 24 )  = (0.1, 0.03, 0.035), (a 31 , a 32 , a 34 ) = (0.3, 0.08, 0.25), (a 41 , a 42 , a 43 ) = (0.01, 0.005,  0.002), (q 1 , q 2 , q 3 , q 4 ) = (0.9, 0.9, 0.8, 0.7), (e 1 , e 2 , e 3 , e 4 ) = (0.005, 0.006, 0.009,  0.015), (d 1 , d 2 , d 3 , d 4 ) = (0.001, 0.002, 0 .005, 0.002) with the initial conditions (x 1 (t 0 ) = 0.01, x 2 (t 0 ) = 0.062, x 3 (t 0 ) = 0.07, x 4 (t 0 ) = 0.1).

Impact of fractional orders in stability analysis-RH stability test 

In this subsection, we analyze the bionomic equilibrium points of the LV system (2) with the harvesting effect applying to all the species. Let c j and p j denote the harvesting cost per unit effect and price per unit of x j , j = 1, 2, 3, 4, respectively. Therefore, Γ is the net revenue function defined by Γ = 4 j=1 ( p j q j x j − c j )e j = 4 j=1 Γ j . Here, Γ j = ( p j q j x j − c j )e j represent the net revenues for respective x j , j = 1, 2, 3, 4. The bionomical equilibrium point (x 1∞ , x 2∞ , x 3∞ , x 4∞ , e 1∞ , e 2∞ , e 3∞ , e 4∞ ) is obtained from the following simultaneous equations:

Here, we discuss the existence of the bionomic equilibrium, which can be done in the following two cases (trivial and non-trivial). If c j > p j q j x j , then the harvesting cost is greater than the revenue of each x j , j = 1, 2, 3, 4. This means harvesting effect is not applied on the species x j , (e j = 0). If c j < p j q j x j , then the harvesting cost is less than the revenue of each x j , j = 1, 2, 3, 4. This means harvesting is applied on the species x j , j = 1, 2, 3, 4.

Case 1:

If the harvesting cost of the all species does not exceed the revenue cost, then the harvesting effect is applied to all species x j , j = 1, 2, 3, 4. Hence, 

Unstable 1

x 2 x * 3 =0.0000, x * 

If the following four feasibility conditions are hold, then e 1∞ > 0, e 2∞ > 0, e 3∞ > 0 and e 4∞ > 0. 

If the harvesting cost of the all species exceeds the revenue cost, then the LV system (2) will be closed. Thus, the above discussion gives the following result.

The LV system (1) has a positive non-trivial bionomic equilibrium point 

Using the Theorem 7, the existence of the non-trivial interior equilibrium point and its feasibility conditions (12) to (15) of the LV system (2) is presented in Table 3 .

Here, we discuss the optimal harvesting policy by using the Pontryagin's maximum principle. Define the cost function J(e 1 , e 2 , e 3 , e 4 ) = t t 0 e −γ t Γ (x 1 , x 2 , x 3 , x 4 , e 1 , e 2 , e 3 , e 4 )dt. Here, γ represents the instantaneous annual discount rate. To maximize the cost function J(e 1 , e 2 , e 3 , e 4 ) subject to the LV system (2) and the control variables e j satisfying 0 ≤ e j ≤ e max j , j = 1, 2, 3, 4. Here, e max j is the feasible upper limit of harvesting effect e j , j = 1, 2, 3, 4. Now, construct the Hamiltonian function H,

Here, λ j (t) = λ j , j = 1, 2, 3, 4 are the adjoint variables. From (16) , for j = 1, 2, 3, 4, we get,

Here, j is called the switching function. The optimal control should be a combination of the extreme control and the singular control as Hamiltonian H is a linear in the control variables e j . Clearly, the optimal control e j that maximizes H satisfies,

Consequently, the optimal harvesting policy satisfies

Here, e * j is called singular control of the optimal control e j , which satisfying 0 < e * j < e max j . According to (17) , we obtain the following equations, for j = 1, 2, 3, 4,

To determine the optimal harvesting effect, introduce the corresponding adjoint equations, dλ j dt = − ∂ H ∂ x j , for j = 1, 2, 3, 4. From Eq. (16), we get, for i, j = 1, 2, 3, 4, dλ j dt = e −γ t p j q j e j + λ j l j r j

With the help of Eqs. (11) , and substituting Eqs. (18) in (19) , and integrating, for i, j = 1, 2, 3, 4,

Here, we omit the integration constants to ensure that λ j e γ t are bounded for j = 1, 2, 3, 4 as t → ∞. Substituting (18) in (20), for i, j = 1, 2, 3, 4, we get the following control parameter e = (e 1 , e 2 , e 3 , e 4 ), e j = h j (γ + a j j l j x j ) i = j q i − i> j h i a i j l i x j i = j q i + i< j h i a i j l i x j i = j q i p j q 2

Here, h j = p j q j x j − c j , j = 1, 2, 3, 4. The optimal harvesting effect e = (e 1 , e 2 , e 3 , e 4 ) can be determined by Eqs. (21) .

We have analyzed the optimal harvesting policy by applying Pontryagin's maximum principle. From the economic and biological aspects of renewable resource management, the harvesting effect is used for all four species x j , j = 1, 2, 3, 4. For environment management, planning of harvest strategies, and keeping sustainable growth are needed. In the ecosystem, we have been concerned with the harvesting effect as a control parameter and its optimal point is e = (e 1 , e 2 , e 3 , e 4 ). If the harvesting effect is applied less than the optimal point e = (e 1 , e 2 , e 3 , e 4 ), then all species will coexist. It can control the optimal level and environmental balance. If the harvesting effect is used greater than the optimal point e = (e 1 , e 2 , e 3 , e 4 ), then the ecosystem will go to dangerous destruction.

In this work, the new problem of the conformable fractional-order four species LV system is built with both immigration and harvesting effects. The existence and uniqueness of the system solution, uniform boundedness, persistence, permanence, and extinction of the LV system (2) are obtained. The stability of the system is obtained with the RH stability criterion. This work verified the small changes in fractional-orders, and the immigration effect gives a significant impact on the dynamic results of the LV system (2) . The numerical results are gained for the obtained theoretical results. The optimal harvesting policy of the LV system (2) is obtained at the bionomic equilibrium point by Pontryagin's maximum principle. Nowadays, researchers are concentrating their study on the mathematical modeling of infectious diseases known as the SIR (Susceptible-Infected-Recovered) epidemic model. In the future, it is motivated that the fractional-order SIR infection model as an essential research. Also, still it is open for stabilization of conformable fractional-order time-delay population system in this direction.

Global dynamical properties of Lotka-Volterra systems

An introduction to the fractional calculus and fractional differential equations

Control systems

The parameter set-1 with the fractional orders α 1 = 0.8, α 2 = 0.95, α 3 = 0.85, α 4 = 0.5 does not satisfy both the RH stability test in Theorem 6 (see Table 1 ) and the condition of the positive interior equilibrium point E(x * 1 , x * 2 , x * 3 , x * 4 ). So that the LV system (2) is unstable (see Fig. 1 ). The parameter set-1 with the fractional-orders α 1 = 0.95, α 2 = 0.79, α 3 = 0.9, α 4 = 0.8 satisfies both the RH stability test in Theorem 6 (see Table 1 ) and the condition of the positive interior equilibrium point. So that the LV system (2) is stable (see Fig. 2 ). The parameter set-1 with the integer-orders α 1 = α 2 = α 3 = α 4 = 1 satisfies both the RH stability test in Theorem 6 (see Table 1 ) and the condition of the positive interior equilibrium point. So that the LV system (2) is stable (see Fig. 3 ). From this analysis, the small changes in the fractional orders of the system can stabilize and destabilize the LV system. Hence, the fractional orders having much influence on the dynamical behavior of the LV system (2) (see Table 1 ). The dynamical and phase-space representations of the unstable and stable behaviors of the LV system (2) are graphically represented in Figs. 1, 2, 3.

To analyze the impact of immigration on stability of the LV system (2) Remark 8 From this case study of the immigration effect on each species of the LV system (2), in the following cases 2, 3, 5, 6, and 7, the LV system (2) is unstable. Since the parameter set-2 with the respective various set of immigration parameters (I 2 , I 3 , I 4 ) does not satisfy conditions of both the RH stability test in Theorem 6 (see Table 2 ) and positive interior equilibrium point. In cases 1 and 4, the LV system (2) is stable (see Table 3 ). Since the parameter set-2 with the respective set of immigration parameters (I 2 , I 3 , I 4 ) satisfying both the RH stability test in Theorem 6 (see Table 2 ) and the condition of the positive interior equilibrium point. From this case study on the immigration effect, we can conclude that the small changes in the immigration effect have a significant impact on the dynamic behavior of the LV system (2) . The dynamical and phase-space representations of the stability (unstable or stable) behaviors of the LV system (2) are graphically represented in figures from Figs. 4, 5, 6, 7, 8, 9, 10.