Population Dynamics of Two Mutuality Preys and One Predator with Harvesting of One Prey and Allowing Alternative Food Source to Predator

In this paper, the interactions among three species populations are considered. The system includes two mutuality preys and one predator. The second prey is harvested. While dependent on preys, the predator has an alternative food source also. The three species interaction can be described as a food chain in which two preys help each other but the predator attacks both the preys according to type I and II functional responses respectively. These population interactions are modeled mathematically using ordinary differential equations. It is shown that the solution of the model is both positive and bounded. The equilibrium points of the model are found and they are analyzed to identify a threshold that will guarantee the coexistence of the populations. Positive equilibrium points of the system are identified and their local and global stability analysis is carried out. Numerical simulation study of the model is conducted to support the results of the mathematical analysis. It is pointed out that as long as harvesting rate on the prey population is smaller than its intrinsic growth rate the coexistence of the system can be achieve. The results of the analysis and the discussion of the population dynamics is lucidly presented in the text of the paper.


Introduction
The interaction created among some communities and the environment is known as ecosystems. The community is a collection of various populations in a given area those interact and influence one another for their survival. In an ecosystem, there will be an interrelation among the organisms in terms of their food chain. The food chain of the populations plays an important role in order to keep a balance in the ecosystem. If one of the foods in the chain is missing then the balanced dynamics of the ecosystem will be disturbed and the populations of preys or predators or both in the ecosystem will be affected.
Mathematical model is the best way used to understand the dynamics and behavior of the ecological interactions among the species such as prey and predator populations. The simple model of Prey -predator interactions is proposed by Lotka and Volterra separately but currently the model is known as Lotka -Volterra model.
The study of Prey -predator species interaction plays an important role in understanding the dynamical behavior of the populations in biological ecology. The selected ecological population systems are not always modeled using a linear system but well modeled through a nonlinear system of equations.
A literature survey on the models of two preys and one predator is conducted as it is the topic of the present research and the review is presented here. The first simple mathematical model that represents the interaction of such three species system and determines their dynamical behavior has been considered and analyzed by [6].
However, these investigational results have been modified Harvesting of One Prey and Allowing Alternative Food Source to Predator by the authors [7]. Also, by including the environmental parameters and on applying algebraic manipulations the criterion for stability of the population sizes in ecology is created in [2]. Further, a functional response mechanism has been introduced. This mechanism expresses the prey density as a function of prey or predator or both. Here, the prey density is defined as the number of preys available per one predator per unit time [11].
Over years various features like coexistence, stability, persistence and extinction as well as predator switching between prey population species have attracted the attention of the researchers. Including the foregoing features and applying Holling type-I functional response a two preysone predator system was addressed by [5].
The existence and global stability of positive periodic solutions of these three species models with the inclusion of Holling type II functional response are studied. Based on the environmental parameters, Holling type II functional response and without the team help the conditions for local as well as global stability are established by [14].
The complexities in the dynamical behavior of two preys and one predator system following Holling type II functional response with impulsive effect is discussed in [3,4,[8][9][10].
It is quite obvious to believe that living in a team or in a group is very important for any population species. The major advantages of living are as follows: (i) Searching for food in a group will be more efficient than doing it alone (ii) Higher probability for predators to attack preys and (iii) Preys defend themselves from their enemies by team.
On the other hand, the local and global stability analysis of these models with the inclusion of group help and Holling type I functional response is carried out in [1]. Here, the three species food web model consisting of two teams of preys interacting with one team of predator is developed based on the assumptions (i) during predation the members of both groups of preys will help each other (ii) the predators are benefitted from predation according to Holling type II functional response (iii) Prey teams help each other for the prevention of predator's attack which will increase the growth rate of prey species [13].
An interaction among two mutualistic preys and predator population has been considered by [12]. Here, the populations were interacts into two areas: free area and refuge area. In free area, only the second prey and predator population species exist and interact while in a refuge area, only the first prey population species exists. In a refuge area, the predator population species cannot enter and attack the prey species. However, in a refuge area the two preys can interact and help each other. In addition, in this model, proportional harvesting function and functional responses are considered among these populations interactions. Based on the unique and positive equilibrium points, local and global stability of the system has been determined analytically and numerically.
In the present paper, the authors have considered interactions among three populations: two preys and one predator. The two prey populations are assumed to help each other. Further, the predator population can attack both the prey populations and gain benefit according to the types I and II functional responses. Additionally, application of proportional harvesting function to one of the prey species is also considered.

Assumptions of the Model
In order to construct a mathematical model, the following assumptions have been made in the present study: (i) All the three species populations grow following logistic function. (ii) There is a mutuality interaction between the first and second prey species. (iii) The first prey population helps the second according to type I functional response while the second helps the first according to type II functional response. (iv) The predator population attacks the first prey population and gets benefit according to type I functional response. (v) The predator population attacks the second prey population and gets benefit according to type II functional responses. (vi) The predator population has alternative food sources.
That is, it doesn't fully depend on any of the prey populations for survival. (vii) The second prey population is harvested proportional to its density. The interactions among the three population species of the present model are visualized through a flow diagram given in Figure 1. In Figure 1 , and represent two preys and a predator respectively. The positive or negative signs indicate that the populations are benefited or affected.

Variables and Parameters
The variables used in mathematically formulating the model assumptions are described here under: (i) is the density of the first prey population at time (ii) is the density of the second prey population at is the density of the predator population at time It is true that , and are time -dependent functions. However, for the simplicity these are represented as , and in the text of this paper. Similarly, the parameters used in mathematically formulating the model assumptions are described here under: (i) The intrinsic growth rates of the first and second preys and the predator populations respectively are , for = 1, 2, 3. (ii) The positive impact on the first prey by the second and vice versa are and respectively. (iii) The helping rate by the second prey to the first is ℎ. (iv) The parameters for = 1, 2, 3 are the carrying capacities of three population species respectively.
(v) The parameters and are the catchability and the effort applied on the second prey population i.e., is the harvesting rate of the second prey population. (vi) The parameter is the negative impact on the second prey population due to the predator population.
(vii) The parameter is the positive impact to the predator population due to the second prey population. (viii)The parameter is the attacking rate of the predator population on the first prey population. (ix) The parameter is the benefiting rate of the predator population from the first prey population.

Mathematical Formulation of the Model
The mathematical model of three population species interaction with the inclusion of functional responses and harvesting function is formulated and is expressed as a system of three ordinary differential equations as given in (1) to (3).

Normalization of the Model Equations
It is well accepted that the normalization of model equations is a very important aspect in dealing with biological models whenever the number of parameters used in the original model is very large. Therefore, normalization of a model reduces the number of parameters. Further, it helps to focus on useful parameters and makes the mathematical analyses simpler.
In order to obtain the normalized or reduced form of the system of equations (1) to (3) the transformation variables employed are = 6, = 7, = 8, = 1 ⁄ : . Thus, the original system of equations (1) to (3) can be expressed in the normalized form as (4) to (6).
In the scaled equations 4 − 6 the notations used

Positivity of the Solution
Here, the constant 6 H denotes the density of the initial population of the first prey at time : = 0 and is a positive quantity by assumption. Furthermore, for any exponent value the exponential function is always positive. Therefore, the solution of 6 : is a positive quantity for all :. Positivity of 7 : : Similarly, to show positivity of 7 : it is appropriate to begin with equation (5) given by K7 K: ⁄ = 7LE − E 7 + D 6 − E − ME 8 1 + E7 ⁄ NO . Equivalently, it can be rearranged as K7 7 ⁄ = LE − E 7 + D 6 − E − ME 8 1 + E7 ⁄ NO K: . In this differential equation, the variables 6, 7 and 8 are the functions of :. Thus, by applying integration and performing some algebraic manipulations, the solution of this differential equation can be obtained as Here, the constant 7 H denotes the density of the initial population of the second prey at time : = 0 and is a positive quantity by assumption. Furthermore, for any exponent value the exponential function is always positive. Therefore, the solution of 7 : is a positive quantity for all :.
Positivity of 8 : : Similarly, to show positivity of 8 : it is appropriate to begin with equation (6) given by K8 K: ⁄ = 8LF − F 8 + G 6 + ME 7 1 + E7 ⁄ NO . Equivalently, it can be rearranged as K8 8 ⁄ = LF − F 8 + G 6 + ME 7 1 + E7 ⁄ NO K:. In this differential equation, the variables 6 , 7 and 8 are the functions of : . Thus, by applying integration and performing some algebraic manipulations, the solution of this differential equation can be obtained as 8 : Here, the constant 8 H denotes Harvesting of One Prey and Allowing Alternative Food Source to Predator the density of the initial population of the predator at time : = 0 and is a positive quantity by assumption. Furthermore, for any exponent value the exponential function is always positive. Therefore, the solution of 8 : is a positive quantity for all :. Thus, from the above three verifications it is very clear that the solutions 6 : , 7 : and 8 : of the model equations (4) -(6) are positive quantities for all : ≥ 0 in the region.

Existence of Equilibrium Points and Their Analysis
On imposing the conditions K6 K: ⁄ = K7 K: ⁄ = K8 K: ⁄ = 0 on the system of model equations (4) to (6)
Locally unstable at \ = 0, 0, 0 : The community matrix p at this equilibrium point reduces to the form as After solving the character equation |p − |}| = 0 and on applying some algebraic manipulations the eigenvalues of p at \ are obtained as | = 1, | = E −E and | = F . Here, it is clear that all eigenvalues are positive since E −E > 0, F > 0 indicating that the system at this equilibrium point is unstable.
Locally unstable at \ = 1, 0, 0 : At this equilibrium point, only the first prey population exists and the other two species are absent. The community matrix at this equilibrium point takes the form as On solving the characteristic equation |p − |}| = 0 of this matrix, the eigenvalues are obtained as | = −1, | = E + D − E and | = F + G . Recall that all the parameters involved here E , D , E , F , G are positive quantities. Also, E > E since the growth rate of the second prey is greater than its harvesting rate. Hence, it is very clear that both the second | and the third eigenvalues | are positives showing that the system at this equilibrium point \ is saddle point and thus is unstable in general.
Locally unstable at \ = 0, 1, 0 : At this equilibrium point, only the second prey population exists and other two species are absent. The community matrix p at this equilibrium point can be takes the form as On solving the characteristic equation |p − |}| = 0 eigenvalues of the matrix p are obtained as | = 1 + D , | = − E + E , | = LF + ME 1 + E ⁄ NO. The parameters involved here D , E E , F , E , E are all positive quantities. Hence, it is easy to observe that both the first and third eigenvalues are positives. Thus, the system at this equilibrium point \ is unstable.
Conditionally Locally stable at \ g = 0, 0, 1 : At this equilibrium point, only the predator population exists and other two preys are absent. The community matrix at this equilibrium point reduces to the form as Therefore, the model is locally asymptotically stable at the equilibrium point \ g as long as the two conditions (i) E > E − E and (ii) G > 1 are satisfied. Otherwise, it is unstable.
Locally unstable at \ h = 6 h , 7 h , 0 : This is a predator free equilibrium point. The community matrix at this equilibrium point takes the form as Here it can be observed that all the three eigenvalues are negative quantities if the following three conditions are satisfied by the parameters: Thus, it can be concluded that the equilibrium point \ i is stable conditionally.
Conditionally locally stable at \ j = 0, 7 j , 8 j : it is the first prey free equilibrium point. The community matrix p at this equilibrium point takes the form as Here in p, the matrix elements represent the expressions as = 1 + D 7 j − G 8 j , = D 7 j , = E − 2E 7 j − E − ME 8 j 1 + E7 j ⁄ N, = ME 7 j 1 + E7 j ⁄ N, = G 8 j , = ME 8 j 1 + E7 j ⁄ N, = F − 2F 8 j + ME 7 j 1 + E7 j ⁄ N. Here it can be observed that all the three eigenvalues are negative quantities if the following three conditions are satisfied by the parameters: (i) < 0, (ii) U• − + 4 X < + and (iii) + − • − + 4 < 0. Thus, it can be concluded that the equilibrium point \ j is stable conditionally.
Conditionally locally stable at \ k = 6 k , 7 k , 8 k : At this equilibrium point all the three populations of the model co-exist. The community matrix at \ k reduces to the form as Here in p k , the matrix elements represent some expressions as K = 1 − 26 k + MD 7 k 1 + E 6 k x ⁄ N − G ‚ ƒ , K = MD 6 k 1 + E 6 k ⁄ N, K = −G 6 k , K = D 7 k , K = E − 2E 7 k + D 6 k − E − ME 8 k 1 + E7 k ⁄ N, K = ME 7 k 1 + E7 k ⁄ N, K = G 8 k , K = ME 8 k 1 + E7 k ⁄ N, K = F − 2F 8 k + G 6 k + ME 7 k 1 + E7 k ⁄ N. The local stability of this kind of system can be obtained using Routh -Hurwitz criteria. The objective of this criterion is to find the roots of the characteristics equation without actually solving it. According to this criterion, the characteristics equation of the community matrix at the given equilibrium point |p k − |}| = 0 can be written as | + "| + …| + † = 0 Here in (8), " = − K + K + K … = − K K − K K − K K − K K + K K + K K † = − K K K − K K K − K K K + K K K + K K K − K K K Recall that the system is stable if the roots lie on the left half plane. In other word, from equation (8), the system is locally stable if " > 0, † > 0 and "… > †.

Global Stability Analysis
In section 1.8, the community matrix at each equilibrium point is found and the local behaviors of the system at each equilibrium point were determined. However, in this section the global stability of the system is considered for the identified equilibrium points. In order to obtain the global stability, the Lyapunov quadratic function is introduced. Here, before determine the global stability of the system, lets write the system equation as a matrix as ‡ ′ =ˆ ‡ Here in (9), the notation ‡ ′ = 6 ‰ , 7, ′ 8′ ‹ ,ˆ= m, n, ℎ , ‡ = 6, 7, 8 ‹ . In order to change the system of equation into linear form, differentiate the right hand side of the equation of (9) with respect to ‡. Thus, the linearized form of the equation can be represented as ‡ ′ = p ‡ Here, p represents the linearized form of the original equation ˆ ‡ which is also called the community matrix. Therefore, the global stability of the system can be determined from the community matrix at the origin.
In this paper the system of equation is globally asymptotical stable at the equilibrium point \ g , \ i and \ j . But, in this section, the global stability at the coexistence equilibrium point \ k is shown here as the following theorem.
Here, Q is any symmetric positive matrix and without loss of generality it can be selected as Q = I . Thus, on solving Lyapunov equation and after performing some algebraic operations the entries of matrix \ are obtained as: Here, it is straight forward to verify that the six elements \ , \ , \ , \ , \ and \ are all positive quantities.
Recall that a matrix \ is said to be a real symmetric positive definite matrix if the determinant of each of the minor of the matrix p are positive. That is if Hence, the eigenvalues of matrix \ also positive and thus the matrix \ is a positive definite.
Recall that an equilibrium point is said to be globally asymptotically stable if the Lyapunov function • ‡ satisfies the following three conditions on the entire state space: (i) • ‡ is positive definite (ii) its time derivative • ′ ‡ is negative and (iii) |• ‡ | → ∞ as || ‡|| → ∞ Thus, the following result: It is already shown that the Lyapunov function • ‡ satisfies all the cited three conditions in case of the interior equilibrium point. Hence, \ k is globally asymptotically stable based on the above two conditions. The following table shows the general behavior of the given system at each equilibrium points.

Result and Discussions
In this model, it is observed different important equilibrium points at which the system is locally as well as globally stable. The stability of the system is determined based on the conditions at which the equilibrium points of the model equations are exists. For example, at equilibrium point \ g only the predator species can survive. That is, since there is an alternative food source for the predator species, it can survive for long time. This result has been shown in figure 2 of the above graph.
In addition, at the equilibrium point \ i , only the first prey and the predator species can live together. Here, as it is observed from the figure 3 of the numerical simulation, the competition is observed between the first prey species and the predator. In fact, the predator species can eat the first prey species but, since the predator species is not fully depend on the prey species they get a good chance to survive for long time without extinct. The behavior of the system at \ j is similar with that of \ i , the difference here is in \ j the predator species interact with the second prey species. The dynamics of these species have been observed in figure 4 of the above graph.
Finally, it is seen that the three species can live together with some conditions described in this paper. Thus, the alternative food source for the predator species helps the preys species to survive for long period of time without extinct.

Conclusions
In this study, the dynamical behavior of three population species is considered. Of the three species two are preys and one is a predator. In this study, different prey dependent functional responses viz., type I and II, are considered. Additionally, the density dependent harvesting function has been imposed on the second prey. At the positive equilibrium points of the model equations the local and global stability analyses of the system have been conducted using Lyapunov function. The positive equilibrium points \ g , \ i , \ j and \ k are shown to be globally asymptotically stable. The behavior or the model populations at these equilibrium points have been demonstrated both analytically and numerically.