On a Three-tier Ecological Food Chain Model with Deterministic and Random Harvesting: a Mathematical Study

In the present study, we consider a nutrient-autotroph-herbivore ecosystem model where the herbivore species is assumed to have a commercial value. We use a Holling type-II harvest function to model density dependent herbivore harvesting. Stability criteria of the resulting model is investigated both from analytical and numerical viewpoints. The investigation revealed the existence of a number of threshold values of the harvest rate that have a remarkable influence on the system dynamics. Next we incorporate a noise term in the parameter representing harvest rate to model the phenomenon of poaching as random harvesting. The stochastic model is analyzed for exponential mean square stability and the resulting criteria in terms of harvest related parameters obtained. These parameter thresholds could be utilized to develop effective harvesting strategies and wildlife management policies which take into account the overall survival of the ecological populations.

On a three-tier ecological food chain model with deterministic and random harvesting: A mathematical study B. Mukhopadhyay 1 , R. Bhattacharyya 2

Introduction
The interrelationship between different trophic levels in an ecological food-chain has interested theoretical and mathematical ecologists for a long time.Understanding the dynamical aspects of these relationship through the formulation and analysis of mathematical models has occupied a significant volume in mathematical ecology.The basic trophic levels in an ecological food chain comprise of (i) primary producers, (ii) primary consumers.The primary producers, also known as autotroph, are capable of producing their organic food requirement from inorganic materials either be photosynthesis or by inorganic oxidation.Herbivores, that feed on these autotroph, are the primary c Vilnius University, 2011 consumers.Ecological food chains basically falls into two categories: (i) the terrestrial plant-herbivore system [1][2][3][4][5]; (ii) aquatic phytoplankton-zooplankton system [6,7].An important and well known aspect in a natural ecosystem is the regeneration of nutrient due to decomposition of dead biotic elements.The effect of nutrient recycling on food chain dynamics has been studied extensively for closed and open ecosystems [8,9].Within the terrestrial environment decomposers like bacteria, fungi etc. break down the organic matter of dead organisms and wastes into inorganic substances.This generates a kind of positive feedback loop in the system that causes a transfer of energy from the higher trophic level into the food resource.In this way, decomposers play a key role in maintaining these resource-based food chains.
Over exploitation of ecological resources to meet growing needs of mankind has been a topic of much concern for ecologists, bio-economists and natural resource managers for some time now.Ever since primitive humans began hunting or fishing thousands of years ago, there has been a need to know how killing a certain number of individuals of a population will affect the population as a whole.The study of population dynamics with harvesting is a subject of mathematical bioeconomics and is mainly concerned with the optimal management of renewable resources [10].Harvesting is commonly practised in fisheries, forestry and wild life management.It has a considerable effect on the dynamical evolution of the harvested species, the severity of which depends on the harvesting strategy that can result from rapid depletion to complete preservation of the concerned population.
The most basic form of harvesting used in mathematical models is constant rate harvesting, that is, Harvesting is modelled by the term h > 0 which indicates that species are harvested at a constant rate h independent of the existing population concentration.An improvement of the above model can be [10], where the catch per unit effort is proportional to the stock P .Here E is termed as the harvesting effort and q the catchability constant.A further refinement is which was first proposed by Holling [11].The effect of constant rate harvesting on preypredator dynamics has been investigated by many authors [12][13][14][15].These investigations revealed very rich and interesting dynamics such as stability of equilibria, existence of Hopf-bifurcation, limit cycles, homoclinic loops, Bogdanov-Takens bifurcation and even catastrophe.Martin and Ruan [16] studied the combined effect of harvesting and time delay on generalized Gause-type prey-predator model and the Wangersky-Cunningham model.Xiao and Jennings [17] performed a detailed bifurcation study of a two species ratio dependent prey-predator model with constant rate harvesting and compared the results with the model with no harvesting.Milner et al. [18] studied temporal and spatial aspects of red deer harvesting in some European countries and explored the biological as well as cultural factors associated with such harvesting.Effects of stochasticity on harvesting in various contexts has also been studied [19,20].
In a real world situation, though there are examples of planned harvesting following some scientific strategy, the majority of harvest phenomena occurs in an unorganized fashion.In case of fishery resources, private ownership is rare except of sedentary species, for example, shell fish, that are readily protected from poachers.In many countries, though off-shore fishing is done in a systematic manner by coastal fishing communities, most of marine fisheries are exploited by competing fishermen and fishing companies [10].In the terrestrial scenario, the situation is much more unpredictable as the herbivore mammals, that have high commercial value, fall prey to poachers in a highly random fashion.Thus, while formulating realistic mathematical models of harvesting, uncertainty in harvesting should be of considerable significance.
In the present research, we take up the task of investigating the implications of herbivore harvesting in an autotroph-herbivore food chain model in the presence of decomposers.As a model system, we take a nutrient-autotroph-herbivore ecosystem foodchain model with nutrient recycling already studied by us in the absence of harvesting [21].We add herbivore harvesting into that model system.We first study the case of deterministic harvesting with a Holling type-II harvest function in an attempt to find out how harvesting affects the system dynamics.Then we add a noise term in the parameter modeling harvest rate and study the stability aspects of the resulting stochastic model.Numerical experiments are performed to complement analytical findings.

Model description
We describe the nutrient-autotroph-herbivore harvesting model by the following three differential equations N (t) stands for the amount of nutrient present within the soil.A(t) denotes the autotroph biomass and H(t) the number of herbivore present at any instant of time t.All these state variable are assumed to be in nutrient equivalent units.We take the initial conditions as We consider external nutrient input to the system at a constant rate N 0 .It is assumed that input of external nutrient is dependent on the amount of nutrient present in the system.
The nutrient uptake rate per unit biomass of autotroph per unit time is α.Nutrient involved in the system also undergo biomass loss due to leaching at a rate a. Autotroph growth, at a rate α 1 , is due to internal positive feedback (photosynthesis in case of phototrophic and inorganic oxidation in case of chemoautotrophic) and they loose biomass at a rate c 1 due to grazing, litter fall etc.Herbivore grazing is modeled using a Michaelis-Menten Holling type-II functional response with β as the grazing rate and q the resulting autotrophherbivore conversion rate.We further assume that herbivores suffer losses due to natural death at a rate b and also due to harvesting.γ 1 and γ 2 are the nutrient regeneration rates from dead autotroph and herbivore population respectively.Harvesting is represented by a Holling type-II function with h as the harvest rate (also known as catchability constant) and E is a measure of the effort required to harvest the herbivores.This type of harvest function implies that when fewer herbivores are available, it is harder to find them and so the daily catch drops.On the other hand when there are sufficiently many herbivores, then lim H→∞ hH E+H = h so that the harvesting level is close to h, the catchability constant.

Boundedness and stability criteria
We start by establishing the biological validity of the model system.
Theorem 1.All non-negative trajectories of (1) which start in R + 3 are uniformly bounded.
The proof is straightforward and is deferred in the appendix.System (1) possesses the following steady states: 1.The trivial state E T ≡ (0, 0, 0).

The axial state E
where Γ = (k+x)h (qβ−b)x−kb and The existence criteria for the boundary equilibrium point A study of the existence criteria of E * reveals that Thus, the harvest rate h is an important parameter in controlling the existence criteria of the interior steady state.
Let us now study the stability criteria around E * .The community matrix around E * is given by The characteristic equation corresponding to J E * is where If a 12 > 0, then simple algebra reveals that B > 0 which in turn implies that D < 0. Consequently, the system can never be locally asymptotically stable in this case.On the contrary, if a 12 < 0, the system may become locally asymptotically stable around E * .Thus depending upon system parameters, the system may exhibit stable or unstable behavior in this case.

Stochastic extension
In this section, we study the effect of random harvesting on the food chain model.To include this effect, we perturb the system stochastically by adding a noise term in the parameter modeling the harvest rate (h).The modified model takes the form where η(t) is a stochastic process that represents the noise term.The precise properties of the noise are not known.However, to gain insight into the effect of noise on the system, we assume the noise to be Gaussian distributed white noise with a zero mean and a delta function auto correlation so that [22] where • denotes expectations.We consider a solution of (9) in the form Substitution of ( 11) into ( 9) yields Again substituting y i (t) = y * i +ξ i (t) and linearizing about (N * , A * , H * ) = (e y * 1 , e y * 2 , e y * 3 ) we obtain where ξ = (ξ 1 , ξ 2 , ξ 3 ) are stochastic perturbations around (y * 1 , y * 2 , y * 3 ).Using the definition of Gaussian white noise η(t) as the derivative of the Wiener process W(t), we write (13) as on t ≥ 0 with initial value ξ(0) = ξ 0 ; where We consider a solution of ( 14) in the form [23] Next we suppose that M has eigenvalues with negative real parts.In that case, we can find a pair of positive constants β 1 and λ 1 such that for t ≥ 0. We further assume that 2 ≤ β 2 e −λ2t for t ≥ 0 for a pair of positive constants β 2 and λ 2 .Now, This implies where λ 1 ∧ λ 2 denotes the minimum of λ 1 and λ 2 .Let, 0 < ε < (λ1∧λ2) 2 be arbitrary.Setting we find that E(| ξ(t)| 2 ) ≤ Θe −(t){(λ1∧λ2)−ε} .Therefore, the system will be exponentially stable in mean square.Now, which is possible when P 2 is small.However, small value of P 2 implies that H * (E+H * ) 2 is small which is possible when h, the harvest rate, is large.
As in our analytical study, we are interested in the interior steady state.This interior state is evaluated as * = (0.8216, 3.55, 2.118).From our numerical study we have obtained three threshold values for the harvest rate h, namely, h 1 , h 2 and h 3 .When h < h 1 (= 0.05), irregular oscillatory behavior in population concentration is observed (Fig. 1(a)).On the other hand for h 1 < h < h 2 (= 0.115), populations evolve in a regular periodic coexistent fashion (Fig. 1(b)).When h is increased further so that h 2 < h < h 3 (= 0.13) stable coexistence occurs (Fig. 1(c)).But as h crosses h 3 , a sudden extinction of the herbivore follows (Fig. 1(d)).Next we simulate the model with noise-induced harvesting using the Euler-Maruyama discretization scheme.We use the same set of parameter values as in the deterministic case with D = 0.0816, the intensity of Gaussian white noise.We first simulate the stochastic model with h = 0.13 for which the deterministic model exhibit stable coexistence.The resulting plot revealed stochastic fluctuation in herbivore time density (Fig. 2(a)).We have also experimented with higher values of the harvest rate and found that for h = 0.14 and h = 0.15 the herbivore exist in an oscillatory fashion (Figs.2(b)-2(c)).However, for h = 0.16, the population undergoes a sharp decline in concentration indicating a possible extinction (Fig. 2(d)).

Concluding remarks and discussion
Harvesting is one of the complex set of factors that influence the vital rates of harvested ecological populations; hence, it would not be realistic to asses the ecological effects of harvesting in isolation [24].Consequently, advancement in the research on harvesting requires a transition from single species models towards a more ecosystem based approach.In line with the above message, we, in the present research, investigated a bio-economic model that deals with a nutrient-autotroph-herbivore food-chain with harvesting of the commercially viable herbivore.The salient features of the model are: (i) density-dependent nutrient input; (ii) nutrient recycling both from autotroph and herbivore; (iii) harvesting of herbivore and protection of autotroph; (iv) use of a density dependent harvesting function.The analytical and numerical study revealed that the harvest rate h played a major role in regulating the system dynamics.We have deduced three critical harvest rates, namely, h 1 , h 2 , h 3 such that when the harvest value crosses these thresholds, the system dynamics changes from chaotic to periodic, from periodic to stable coexistence, and from stable coexistence to herbivore extinction.Interestingly, in one of our previous studies involving the same food-chain model without harvesting [21], the herbivore mortality rate was seen to be the controlling parameter; and it was observed that for small mortality of herbivore, species coexistence occurs, whereas for large mortality, the herbivore will become extinct.
In the recent past, the crime of poaching on animals having considerable commercial value has increased beyond proportion.As some of these animals are in the endangered category, such events are creating havoc on the ecosystem as a whole.In the Indian context, the one-horned rhino, the musk deer, the black buck deer, the tusker are to name a few.Due to the presence of enforcing laws and time varying security measures in various sanctuaries, poaching rate fluctuates in a highly random fashion.We have modeled this rate by adding a noise term to the harvesting rate.We have derived the criteria for exponential mean square stability for the stochastic model and showed that in this case also the harvest rate h controls stability aspects.Numerical simulation of the stochastic model revealed another interesting result -the harvest rates that ensure stable coexistence (h = 0.13) and herbivore extinction (h = 0.14) in the deterministic case, induce instability in the stochastic context.However, a further increase in the said rate (h = 0.15-0.16)has a stabilizing effect and a subsequent extinction tendency on the herbivore population.Thus, the harvest rate that causes stabilization and subsequent extinction in the stochastic context is higher than that in the deterministic one.
To summarize, our analytical and numerical investigations revealed the existence of a number of important thresholds for the harvest rate in relation to the deterministic and the stochastic model which could be utilized by field ecologists in formulating optimal harvesting strategies and in developing wildlife management policies that take into account the overall survival of endangered species.

Appendix: Proof of Theorem 1
We take W = N + A + H. Then is the region of attraction and consequently all non-negative trajectories of (1) that initiate in R + 3 will be uniformly bounded, that is, every non-negative trajectory which is defined on [0, ∞) is bounded on an arbitrary compact subinterval of it.