Analysis of Prey – Predator System with Prey Population Experiencing Critical Depensation Growth Function

In this paper we have presented a pair of coupled differential equations to represent a prey – predator system. It is assumed that the growth of the prey population follows critical depensation function and that of the predator population is negative in absence of the prey population. The critical depensation function is special since the growth rate is negative initially but positive later on. This function is stable both at the origin and at the carrying capacity while unstable at the critical mass quantity. The maximum and minimum rates of the critical depensation model are verified. It can be interpreted here that the prey represent fish and the predator represent a kind of birds that mostly feeding on fish to live. We showed the solution of the model is positive and bounded. The mathematical model of the system consisting of 7 parameters is constructed and shown that the non – dimensionalization decreases the number of model parameters to 4. The deterministic behavior of the model around feasible equilibrium points and criteria of the interior positive equilibrium points and their stability are explained. The trivial equilibrium point is always stable while the two axial equilibrium points and the lone interior equilibrium point are either stable or unstable depending on the conditions imposed on the parameters. The criterion for the existence of the limit cycle and the region of existence of interior equilibrium point are identified. Global stability of interior equilibrium points is also studied. For the interior equilibrium point of the model (i) the region of existence is identified (ii) Dulac’s criteria is applied to find the limit cycle and (iii) Lyaponov function is used to analyze the global stability. Simulation study of the model is conducted in support of the analytical analysis. To solidify the analytical results numerical simulations are provided for hypothetical set of parametric values.


Introduction
Mathematical modeling is a key tool that has been considered and used for understanding both growth and dynamics of a population [1 -3]. Mathematical models are widely used for understanding the present size and growth rate of a population and to estimate both past and future population sizes. Mathematical modeling plays an important role in managing any industry, for example fishing industry or a fishery. Also in case of a park where two or more types of populations or species live, mathematical modeling can be applied so that a balance among the populations can be maintained. Further, population models are used to determine maximum harvested quantities in case of agriculturists [4].
The relationship between the population sizes of species in a system is highly non-linear. However, mathematical modeling equations are being designed and used to understand that relationship. Research in the area of theoretical ecology was initiated in 1925by Loteka and Volterra [5]. The Lotka-Volterra predator-prey model is one of the earliest prey predator models which is based on basic mathematical logic [6]. In this paper we have studied the prey -predator model replacing the classical model's exponential growth of the prey population by critical depensation growth with carrying capacity and critical mass quantity . In the present model the predator has no alternative food except fish as that of the generalized preypredator model [7]. The predicted stability of the present mathematical model is well supported by the simulation study also.

Lotka -Volterra Predator -Prey Models
The classical model describing the predator-prey problem was originally developed in the 1920s by Vito -Volterra, and since then a wide and related studies have been conducted. The Lotka -Volterra predator-prey equations are first-order and nonlinear differential equations defined as [1,2,3,4,6,7,9].
In the classical model as described in (1) the variable represents the size or density of a prey population while represents that of a predator population. The quantities , , and are all positive parameters.
The equation (1a) describes the population dynamics of prey in the prey -predator environment. The growth rate of prey ⁄ is proportional to its own population and the declination rate is proportional to the interaction between prey and predator . Here is intrinsic growth rate of the prey population while is the predation rate of prey. The equation (1b) describes the population dynamics of predator in the prey -predator environment. The death rate of predator ⁄ is proportional to its own population and the growth rate is proportional to the interaction between prey and predator . Here is death rate of the predator population while is the growth rate of predator due to interaction with prey.
Equilibrium points and stability. Population equilibrium points occur for the model (1) when growth rates of population are vanishing. It is well known that the equilibrium points for the model (1) are 0 , 0 and ⁄ , ⁄ . Using non-linear stability analysis we conclude that the trivial equilibrium point is unstable while the interior equilibrium point is a center which is stable.
The modified Lotka -VolterraPredator -Prey Models. The logistic inclusion predator prey modelis more realistic than the the previous model and it is written as Here is intrinsic growth rate, is the rate at which the predator and prey meet, is death rate, is the rate at which the predator population grows and is the carrying capacity of the prey population.
Equilibrium points and stability Population equilibrium points occur for the model (2) when growth rates of population are vanishing. It is well known that the equilibrium points for the model (

The Critical Depensation Model
The population growth models can be classified into three types namely Compensation, Depensation and Critical depensation [13]. Compensation growth is a growth type where population declination is compensated by increased growth rate and it is written in the form ⁄ = 1 − ⁄ . Depensation growth is the opposite case to compensation growth model and has equation of the form ⁄ = 1 − ⁄ . In both of these model s denotes the population size at time , is the growth rate, is the carrying capacity and ≠ 1 is any real number different from a unit. But the critical depensation model is extremely in opposite of the depensation model and is briefly explained here in what follows.
As it is clearly described [9] some populations experience reduced rates of survival and reproduction when reduced to very low densities. This reduced per -capita growth rate at low densities is called depensation. The strong Depensation is called critical depensation. Mathematical expression of the strong critical depensation is given by the cubic growth model as, Here in the growth model (3), represents the population size, represents intrinsic growth rate of population, is the carrying capacity of the population in the environment, is the critical mass quantity and ⁄ is the growth rate of the population. In case of the critical depensation model (3) the following observations can be made.
i. Growth rate of the population is negative as long as the population size lies below the critical mass quantity . That is, the condition ⁄ < 0 satisfies in the interval 0 < < . ii. Growth rate of the population is positive as long as the population size lies above the critical mass quantity and below the carrying capacity . That is, the condition ⁄ > 0 satisfies in the interval < < .
iii. Both the critical mass quantity and the carrying capacity are positive quantities such that 0 < < .
iv. The per capita growth rate is always a positive quantity ⁄ ⁄ > 0.
v. The per capita population growth rate ⁄ ⁄ increases with increasing over some range of population size 0, . vi. The critical depensation model (3) has three equilibrium points = 0 , = and = respectively. The equilibrium point = is unstable and the remaining two equilibria points = 0 and = are asymptotically stable. The geometrical representation of these equilibria points is shown in Figure 2. vii. If the initial population size is assumed to start with some value above then the population will grow and reach the carrying capacity over time otherwise it will die out over time. viii. Together with the natural positive population growth if either harvesting or hunting is introduced then the population size can be brought down. As long as the population size is maintained above then the population naturally grows and reaches the carrying capacity . That is, the population size continues growing and reaches the carrying capacity as long as the harvested or hunted population quantity satisfies the following inequality:

"
Harvested or Hunted size 0 < 1" Present population size 0 − " Critical Mass Quantity k 0; ix. If the population size is brought below due to harvesting or hunting then the population size naturally decreases and dies down eventually. That is, the population size continues decreasing and reaches extinction whenever the harvested or hunted population quantity satisfies the following inequality: Implicit solution of the critical Depensation model. An implicit solution of the model (3) is found by applying separation of variables method. The solution is called implicit since the population size variable is not expressed in terms of the dependent variable . Application of separation of variables on (3) immediately gives us L − − ⁄ = L . The integrand on the left hand side can be simplified using partial fraction decomposition method. Let Here M, N and G are unknown constants and these are to be determined. We now determine the values of the coefficients M, N and G as follows. After simplification and comparison of the coefficients we obtain M 1 , N ⁄ and G ⁄ and thus the integral equation takes the form O ln G represents the integral constant.The implicit solution can be further simplified on applying antilogarithmic function and can be expressed as The result (4) is the required general implicit solution of the critical depensation model described by (2). Up on substituting the initial condition at the general solution (3) we obtain.
The result (5) is the required particular implicit solution of the critical depensation model described by (3). In Figure 1 we have shown the growth curve of the critical depensation model for some particular values of the parameters as shown. The curve is plotted for the population size function versus the population growth rate function Z ⁄ . The growth rate is negative in the interval 0 , while it is positive in the interval , . In the interval , the growth rate increases up from 0 to maximum and then falls down to 0. The maximum growth rate i.e., I Z occurs when the population size assumes the value 3 ⁄ and the corresponding maximum growth rate is given by I Similarly, the minimum growth rate i.e., ICB Z occurs when the population size assumes the value 3 ⁄ and the corresponding minimum growth rate is given by Alternatively there is a geometric approach that is useful for verifying local stability of equilibrium points in case of the first order ordinary differential equation (3). The geometric approach involves graphing the direction of flow on a phase line diagram as shown in Figure 2   is unstable (iii) The equilibrium point = is stable since all the solution curves that start with any initial value around are running towards and approaching = as tends to infinity.

The Proposed Prey -Predator Model
We study the Predator -Prey model considering fish population as prey while birds population as predator. The newly proposed predator -prey model can be described mathematically as In the model (6), the parameter represents the carrying capacity of the prey population, is critical mass quantity of the prey populationsuch that 0 < < , and are intrinsic growth rates, e is maximum up take rate, and H is half saturation value of the predator. As we observe in the proposed prey -predator model (6), the predator cannot survive in absence of the prey. The per capita growth rate of the prey is decreased by the quantity eK H + ⁄ in the presence of predator. The per capita growth rate of the predator is enhanced by the quantity > in presence of prey. The model (6) is designed based on the assumptions as described below.
i. The variable represents the size or density of prey population while K represents the size or density of the corresponding predator population. ii. The prey population increases due to the mutual benefit between the individuals of the prey population with the rate + ⁄ . iii. The prey population decreases (a) with an intrinsic growth rate , (b) with multi species intra specific competition rate ⁄ f among prey and (c) with the rate e K H + ⁄ because of the bird predation given by E = e H + ⁄ .
iv. The predator population dies out in the absence of the prey. v. The predator population increases with rate > K because of the presence of prey population .

Non-dimensionalization of the Model
Scaling decreases the number of parameters of the model and it simplifies the model equations. Let us replace the old variables , K and with a new set of variables A , = and g through the transformation equations given by = A, K = e ⁄ = and = 1 ⁄ g . Up on using these transformation equations of the scaled variables in (6), we obtain Here in (7) the dimensionless expression of the model, q = ⁄ , = H ⁄ , o = ⁄ and p = > ⁄ are dimensionless parameters. The term A + A ⁄ is called functional response term. Further, we note that A , = and g are dimensionless variables and either of them does not depend on the maximum uptake rate e.

Boundedness of Solution of the Model
Here in what follows it is to be shown that the functions We further show that for the model (7) together with positive initial conditions there exists a unique solution. For that, it requires that both n A, = and n A, = are Lipschitz functions within the domain r d = A, = : A > 0, = > 0 . We first show that n A, = is a Lipschitiz function as [9]. Note that a function v , K is said to be a Lipschitz function if the condition |v , K − v , K x | ≤ z |K − K x| is satisfied. Note that we have used the fact | u A ⁄ | 1 . Similarly, the function n A, = is also a Lipchitz function. Therefore, the solution of the model (6) with non-negative initial conditions exists and is unique. We now assume and require that the solution of the model (6) (7) is invariant under the region r d . That is, the solution is positive and remains the region r d [9]. Proof From the prey equation of the model (7) . Thus, it can be concluded that 0 ≤ " g ≤ z o 1 − ⁄ . Therefore, the model (7) is bounded or dissipative with the asymptotic boundary z o 1 − ⁄ . Hence, the solution of the model (7) is bounded or dissipative.
In Theorem 2 we observe that (i) Here we have used the notation ℎ = 1 q ⁄ .

The Equilibrium Points
The equilibrium points of the system (6)  Theorem 2 The prey predator model (7) with prey experiences critical depensation growth has a positive interior equilibrium point ‰ f if and only if 0 < q < A * < 1 or equivalently 0 < q < o p ⁄ < 1 holds true where A * = o p ⁄ . More over Positive interior equilibrium point does not exist for o p ⁄ < q and o p ⁄ > 1 .

Local Stability of the Equilibria
The local stability of an equilibrium point is determined by linear stability analysis of which by observing the algebraic sign of the eigenvalues of the Jacobian matrix. An equilibrium point is locally stable if we perturb the initial condition slightly, then the system stay in the neighborhood of that equilibrium point or for asymptotically stability, the system returns to the equilibrium point [11]. The Jacobin matrix for (7) is constructed as Š A, = = T V .
Here  Proof: As mentioned above the stability conditions are g J E f < 0 and "> J E f > 0. g "> J E f 0 ⇔ 0 q o p ⁄ 1. Therefore after simplification we conclude that the theorem is proved.
Positive interior equilibrium point and the phase diagram In Figure 3, we observe that the prey zero growth isocline = 1 A ~ A q ⁄ 1• A is a cubic function and the predator zero growth isocline pA o 0 is a vertical line. The cubic function intersects the A axis at A q and A 1 and has a positive maximum value < < < f q ⁄ where < ~~q √q 3• 3 ⁄ • i.e. = A ⁄ A 0 if A ~~q √q 3• 3 ⁄ • ≡ < and the corresponding value of = is < 1 q ⁄ < < f 1. The intersection of these two isoclines in the A= plane gives the following facts: i. There is no positive interior equilibrium point for o p ⁄ q and o p ⁄ 1 but that exists for 0 q o p ⁄ 1. ii. The interior equilibrium points are unstable when q < o p ⁄ 1 and they are stable when q o p ⁄ < 1. iii. The limit cycle exists at the point <, < 1 q ⁄ < < f 1 . We now verify whether or not there exists a limit cycle for the model (7) using Dulac's criteria. According to Dulac's criteria [11] a limit cycle exists for (7) in the positive quadrant of A= planeif the expression 1 uv ⁄ such that N A, = 0 for A 0 and = 0. In case of the present model (7)

Global Stability of the Interior Equilibrium Point Ÿ
A powerful tool for determining global stability of an equilibrium point is the use of Lyapunov function. The Lyapunov function ¡ A, = for the model (6) can be constructed [5] as ¡ A, = A A * A * FB A A * ⁄ = = * = * FB = = * ⁄ . Now, the total derivative of the Lyapunov function is given by Therefore, the Lyapunov theorem [5] implies that the interior equilibrium point ‰ f A * , = * is globally asymptotically • andthis further simplifies to the inequality = A ⁄ A * = * ⁄ . Here A * , = * is interior equilibrium point.

Simulation Study
A numerical simulation is carried out for various choices of biologically feasible parameter values and for different set of initial conditions. From the phase diagram given in Figure 4, we observe that the interior equilibrium points exist for q o p ⁄ 1 and these interior equilibrium points are unstable if q o p ⁄ < 1 . The trivial equilibrium point is stable while the axial equilibrium points are unstable. From the phase diagram given in Figure 5, we observe that the interior equilibrium points exist for q o p ⁄ 1 and they are stable for q o p ⁄ < 1 . The trivial equilibrium point is stable while the other axial equilibrium points are unstable. From the phase diagram given in Figure 6, we observe that the limit cycle exists for interior equilibrium point q o p ⁄ < 1. The trivial equilibrium point is stable while other axial equilibrium points are unstable.
In Figure 7, the parametric values satisfy conditions of Theorem 2. We observe from the time series plot the prey population increases slowly and it converges to a finite value. On the other hand the predator population size decreases for some time and it converges. The prey vs Predator plot also shows size of the predator decreases and increases as the prey size increases.

Conclusions
The paper studied the prey -predator model replacing the classical model's exponential growth of the prey population by critical depensation growth with carrying capacity and critical mass quantity . The model is meaningful more if the prey population is fish and the corresponding predator population is a bird since birds are best feeders of fish. We showed that the constructed model has a unique and bounded solution. The interior equilibrium point ‰ f A * , = * is positive if the condition 0 q A * o p ⁄ 1 holds true. The maximum and minimum rate values of the critical depensation model are determined in this study. The criteria for feasibility of equilibrium points, the stability criteria of the interior equilibrium point and the condition for the existence and nonexistence of the limit cycle in the positive phase plane are explained. The limit cycle exists at the point <, < 1 q ⁄ < < f 1 where < ~~q √q 3• 3 ⁄ • . The simulation study supports the stability of the system predicted using the proposed mathematical model describing preypredator system and all parameter values satisfy condition of theorem 2.