Global stability of two-species mutualism model with proportional harvesting

Article history: Received 2 March 2017 Received in revised form 28 May 2017 Accepted 7 June 2017 This paper deals with the study on a mathematical model consisting of mutualistic interactions among two-species with proportional harvesting. Harvesting function is introduced to describe the rate of removal of the species. The local stability analysis shows that the unique positive equilibrium point is asymptotically stable when certain conditions are satisfied. Global stability is discussed by constructing Lyapunov function. It has been shown that the unique positive equilibrium point is globally asymptotically stable. Finally, numerical simulations supporting theoretical results are also included.


Introduction
*Mutualism is an interaction between two species where both species benefit from the interaction (Rockwood, 2015). Some common examples of twospecies mutualism include zebra and wildebeest (Fay and Greeff, 2006), yucca moths and yucca plants (Starr et al., 2015) and damselfish and sea anemone (Benz, 2000). In a thorough review of the natural history of mutualisms, Janzen (1985) has argued that most mutualisms can be classified into four types; seed dispersal mutualisms, pollination mutualisms, digestive mutualisms and protection mutualisms. Some example of each type of mutualisms can be seen in Kot (2001). In general, mutualism may be facultative or obligate (Morin, 2011). In the case of facultative mutualism, each species is able to survive in the absence of the other. Zebra and wildebeest, damselfish and sea anemones are examples of facultative mutualism. A mutualism is obligate when one species cannot live without the other species. Yucca moths and yucca plants have a reciprocal obligate relationship.
One of the important issues in a mutualism system is global stability. Some biologists believe that local asymptotically stable equilibrium point is globally asymptotically stable in an ecological model (Cheng et al., 1981). A general method to establish global stability in a mutualism system is by constructing a Lyapunov function. Two different forms of Lyapunov functions have been constructed by Goh (1979) and Reddy et al. (2011) to show that the unique positive equilibrium point is globally asymptotically stable in a same basic model of mutualism.
Recently, Georgescu et al. (2016) have analyzed the dynamics of three models of mutualism and establishing the global stability via the method of Lyapunov functional. The authors found that the use of higher-order self-limiting terms cures the shortcomings of Lotka-Volterra mutualisms, preventing unbounded growth and promoting global stability.
In the context of mutualism interaction, some studies consider the effect of harvesting and time delay. The removal of a number of a species of population from its habitat is known as harvesting. According to Ouncharoen et al. (2012), harvest management is used to control increasing population and to meet the public demands for animal damage control, recreation or commercial harvesting. A delay is a general concept that can represent different phenomena such as the finite gestation period or maturation period of a species. For animals, the delay time may be the time it takes for an egg to develop into a fertile adult (Haberman, 1998). The global stability of the unique positive equilibrium point of mutualism model subject to proportional harvesting has been studied by León (2012). He constructed more than one Lyapunov function to prove the global stability of the same model.
The global stability analysis of mutualism model with time delay has been conducted by some researchers. Hsu and Ho (2006) obtained the conditions for global stability of two-species mutualism model with single delay by constructing a Lyapunov function. The authors illustrate their results by some example. The global stability analysis of mutualism models with two delays or without delays by constructing a Lyapunov function has been studied by León (2015). He found that if all the delayed feedback is positive then the stability properties are independent of the values of the delays. Saito (2002) derived necessary and sufficient conditions for the global stability of mutualism model with time delays in some specific cases using the method of Lyapunov function.
In this paper, we extend the work done by Goh (1979) and Reddy et al. (2011) by constructing a suitable Lyapunov function to prove the global stability of the unique positive equilibrium point of a mutualism model with proportional harvesting and illustrate our results by a numerical example. This model is different to a model appearing in León (2012). In particular, the mathematical analysis is simulated using the data for the population of wildebeest and zebra as reported in Fay and Greeff (2006).

Lyapunov's second method
Lyapunov's second method or direct method uses an energy-like function called the Lyapunov function to study the behavior of dynamical systems analytically (Do and Pan, 2009). This method is referred to as a direct method because no knowledge of the solution of the system of differential equations is required. Lyapunov's second method enables the analysis to be extended beyond a small region near the equilibrium point (global analysis). The basic idea of this technique for verifying stability is to seek an aggregated summarizing function that continually decreases towards a minimum as the system evolves. The following Definition 1 and Theorem 1 can be obtained from Boyce and DiPrima (1992). Theorem 1: If there exists a Lyapunov function, then the equilibrium point is stable. If furthermore, the function is strictly negative for every point then the stability is asymptotic.

Mathematical model
We consider the mutualism model with proportional harvesting as follows where, = and = . The symbols and denote the population size of the first species and the second species respectively, and are the intrinsic growth rates, and are the carrying capacities, ⍺ and β measure the mutualism effect of on and on respectively, ℎ and ℎ are harvesting terms proportional to the respective population size and all the parameters are positive constants. Let 1 = − ℎ and 1 = − ℎ , so model (1) can be written as follow Under the assumptions 1 > 0 and 1 > 0, model (2) (2) if 1 = − ℎ > 0, 1 = − ℎ > 0 and − > 0. In this case, the range of harvesting level for model (2) are 0 < ℎ < and 0 < ℎ < .

Local stability
In this section, the local stability analysis of equilibrium points of model (2) is investigated. In order to determine the stability of each equilibrium point, we need to compute the Jacobian matrix of model (2). The sign of the real part of the eigenvalues of this matrix evaluated at these equilibrium points determines its stability. The Jacobian matrix of model (2) takes the form The Jacobian matrix of model (2) at E 1 is The eigenvalues λ 1 = r 1 , λ 2 = s 1 of J 1 are both positive. Thus, the equilibrium point E 1 of model (2) is unstable. The Jacobian matrix of the model (2) at The eigenvalues of J 2 are λ 1 = −r 1 , λ 2 bs 1 +βr 1 b . Since λ 1 < 0 and λ 2 > 0, the equilibrium point E 2 of model (2) is a saddle point. The Jacobian matrix of model (2) at E 3 is The eigenvalues of J 3 are λ 1 = er 1 +αs 1 e , λ 2 = −s 1 . Since λ 1 > 0 and λ 2 < 0, the equilibrium point E 3 of model (2)  Under the condition be − αβ > 0, then (bx * + ey * ) 2 − 4(bex * y * − αβx * y * ) > 0 and (be − αβ)x * y * > 0. Thus, both eigenvalues have negative real roots. Since λ 1 < 0 and λ 2 < 0, the unique positive equilibrium point E * of model (2) is locally asymptotically stable.

Numerical simulation
Now, we carry out numerical simulation to investigate the dynamic behavior of the system (1) about the equilibrium point * = ( * , * ). Consider the following system of wildebeest and zebra without the presence of lion as proposed in Fay and Greeff (2006) by incorporating proportional harvesting to illustrate the theory discussed previously where and represent the population size of wildebeest and zebra respectively (both measured in thousands). Harvesting refers to the process of cropping the species, i.e. removing the species to a new location in order to avoid overcrowding.
Comparing model (16)  We would like to determine the stability of the model (16) where the equilibrium point is ( * , * ) = (14.3365,18.5905). In Fig. 1, the equilibrium point (14.3365,18.5905) is a stable sink. This means that solution trajectories tend to the equilibrium point without oscillation. The population dynamics of wildebeest and zebra can be analyzed separately with respect to time. From Fig. 2

Conclusion
This study considered a two-species mutualism model with proportional harvesting. A suitable Lyapunov function was successfully constructed in order to prove the global stability of the unique positive equilibrium point when the conditions 1 = − ℎ > 0, 1 = − ℎ > 0 and − > 0 are satisfied. From the numerical analysis, it is observed that both the species persist. As seen from Fig. 1-Fig.  3, both the wildebeest and zebra populations tend toward the equilibrium point (14.3365,18.5905). There is still a tremendous amount of work to be done in this area. For example, it is interesting to prove the global stability of two or more species mutualism models having both delay and harvesting function in a real life situation. Previous researchers have studied models with either delay (Hsu and Ho, 2006;León, 2015;Saito, 2002) or harvesting (León, 2012) but as far as we know, no research has been done for models having both components.
The literature on the stability of the mutualism model is not as much as the literature concerning the stability of the predator-prey and competition models. Therefore, the findings of this study may contribute to a better understanding of the population dynamics of mutualism model particularly in the use of Lyapunov function to determine the global stability when the model incorporates harvesting function. Continuous effort in this area can add to the body of knowledge and benefits the society especially regards to the conservation of species within the ecosystem.