A Modified Holling-Tanner Model in Stochastic Environment

Abstract. Recently, a modified version of the so called Holling-Tanner model is introduced in the ecological literature. A detailed accoun t of the deterministic dynamics of this model is presented. The growth rates of the prey and pr e ator are then perturbed by Gaussian white noises to take into account the effect of flu ctuating environment. The resulting stochastic model is cultured by the technique of s tatistical linearization and criteria for non-equilibrium fluctuation and stability are d rived. Numerical simulations are carried out. The implications of our analytical findings are addressed critically.


Introduction
Over the last century, mathematical biology research opens up a new exciting cornucopia of challenging problems for the mathematicians.On the other hand, for the biologists, mathematical modelling offers another research tool commensurate with new powerful laboratory techniques.Different mathematical techniques have been successfully derived to get an insight into the problems of biology, but there is no denying that it is a strong tendency of the mathematicians to neglect factors which could mar the beauty of the analysis.In widening and deepening the scope of mathematical biology, consideration of some relevant factors are very important.
The most crucial element in predator-prey models is the "functional response" or "trophic function", the function that describes the number of prey consumed per predator per unit time for given quantities of prey x and predator y.Various forms of functional responses have become the focus of considerable attention from time to time in ecological literature.The so called Holling-Tanner model is concerned with the Michaelis-Menten or Holling Type-II functional response of the form p(x) = cx m+x , where c is the maximal predator per capita consumption rate, i.e. the maximum number of prey that can be eaten by a predator in each time unit and m is the half capturing saturation constant, i.e. the number of prey necessary to achieve one-half of the maximum rate c.Such functional responses are labelled "prey-dependent" by Arditi and Ginzburg [1] since it depends on prey density only.It was recognized early that the predator density could have a direct effect on functional response.A number of such "predator-dependent" models have been proposed, the most widely known being those of Hassell and Varley [2], DeAngelis et al. [3], Beddington [4] and Arditi and Akcakaya [5].A predator-dependent functional response, which is a function of the ratio of the prey and predator, is known as a ratiodependent functional response.Arditi and Ginzburg [1] introduced a Michaelis-Menten type ratio-dependent functional response of the form p(x/y) = c(x/y) m+(x/y) = cx my+x , where x, y stand for densities of prey and predator respectively.The positive constants c and m are the capturing rate and the half capturing saturation constant, respectively.Predatorprey models with such ratio-dependent functional response are strongly supported by numerous field and laboratory experiments (Arditi and Ginzburg [1], Arditi and Berryman [6], Arditi et al. [7], Hanski [8], Arditi and Saiah [9], Gutierrez [10], Blaine and DeAngelis [11], Poggiale et al. [12], Bernstein et al. [13], Cosner et al. [14], Arditi et al. [15]).Detailed arguments on the merits of ratio-dependent model in comparison with other models may be found in the works of Berezovskaya et al. [16] and Arditi et al. [17].Various aspects of the Michaelis-Menten type ratio-dependent predator-prey system and the effect of many relevant factors like time-delay, diffusion, environmental stochasticity etc., have been extensively discussed in the literature; most recent works being those of Bandyopadhyay and Chattopadhyay [18], Tang and Zhang [19], Huo and Li [20], Lizana and Marin V. [21], Maiti and Samanta [22], Wang et al. [23], Zeng [24], Zeng and Zhai [25], Maiti et al. [26,27], Ruan et al. [28].
Most of the models proposed and studied in the ecological literature work within the framework of an unvarying, deterministic environment.However, the parameters characterizing real environments all, to greater or lesser degree, exhibit random fluctuation.That is, real environments are uncertain, stochastic.Therefore, in deterministic situation, it is always difficult to predict the future of the system accurately.One reason to this difficulty is that biological systems are subject to apparently random fluctuations.In fact, randomness or stochasticity plays a vital role in the structure and function of biological systems (May [29], Renshaw [30], Nisbet and Gurney [31], Samanta [32]).The environmental factors are time-dependent, randomly varying and should be taken as stochastic.Renshaw [30] mentioned that the most natural phenomena do not follow strictly deterministic laws but rather oscillate randomly about some average so that the deterministic equilibrium is not an absolutely fixed state; instead it is a "fuzzy" value around which the biological system fluctuates.Both demographic stochasticity and the environmental stochasticity play a significant part in the realistic dynamical modelling of ecosystems.A central obstacle in the stochastic modelling of an ecosystem is the lack of mathematical machinery available to analyze non-linear multi-dimensional stochastic process.A quantum leap in the mathematical sophistication of ecological modelling occurred when May [29] introduced stochastic differential equations to investigate li-mits to niche overlap in randomly fluctuating environment.Subsequently, the sensibility of stochastic models in comparison with deterministic models is established by many researchers (see Dimentberg [33,34], Samanta [32,35], Samanta and Maiti [36,37], Bandyopadhyay and Chakrabarti [38], Bandyopadhyay and Chattopadhyay [18], Maiti and Samanta [22,39], Maiti et al. [26] and references cited there in).Also, the effects of environmental fluctuations on the models with ratio-dependent functional responses have now become the focus of considerable attention in the literature (see Bandyopadhyay and Chattopadhyay [18], Maiti and Samanta [39], Mankin et al. [40], and Maiti et al. [26]).
Recently, a modification of the so called Holling-Tanner model by invoking the ratiodependent functional response is suggested by Haque and Li [41].In this paper, a detailed dynamics of the model is presented when the environment is assumed to be deterministic.Our analytical findings are illustrated through computer simulations.Wollkind et al. [42], Collings [43] have used the Holling-Tanner model to study the population interaction between the predacious mite Metaseiulus occidentalis Nesbitt and its spider mite prey Tetranychus mcdanieli McGregor.It has been reported by many researchers that environmental fluctuations play a very significant role on M. occidentalis -T.mcdanieli interaction (e.g., Wollkind and Logan [44], Wollkind et al. [42], Collings [43]).To take into account the effect of fluctuating environment, we formulate the stochastic version of the model by perturbing the growth rates of the prey and predator by Gaussian white noises.The criteria for non-equilibrium fluctuation and stability are derived.A comparative study of the stability behaviour of the model in deterministic and stochastic environment is presented.

The basic deterministic model
The May or Holling-Tanner model for predator-prey interaction is with Here N 1 (T ) and N 2 (T ) denote prey and predator densities, respectively, in time T .It is assumed that in the absence of the predator, the prey population density grows according to a logistic curve with carrying capacity K and with an intrinsic growth rate (or biotic potential) r.The parameter s denotes the intrinsic growth rate of the predator.c is the maximal predator per capita consumption rate and m is the half capturing saturation constant.The predator growth equation is of logistic type with a modification of the conventional one.Here the available resources is not constant, but is equal to nN 1 , which is proportional to prey abundance.The parameter n is the measure of the food quality that the prey provides for conversion into predator births.Several dynamical behaviours of Holling-Tanner model have been studied extensively in literature (see May [29], Tanner [45], Wollkind et al. [42], Murray [46], Hsu and Hwang [47], Collings [43,48], Sáez and González-Olivares [49], Braza [50]).
It is already mentioned that Haque and Li [41] have introduced a modified version of the above Holling-Tanner model by replacing the Holling type-II prey-dependent functional response with a ratio-dependent one.Although not mentioned clearly, the main reason for such modification would be perhaps the widespread popularity of the ratiodependent functional response in comparison of its prey-dependent counterpart.Further, it seems that Haque and Li [41] might have expected some different result on "paradox of enrichment" (which they have successfully obtained also).However, as a starting point of our study, we take their modified model described under the framework of the following set of ordinary differential equations: with The model we have just specified has six parameters, which makes analysis difficult.To reduce the number of parameters and to determine which combinations of parameters control the behaviour of the system, we nondimensionalize the system (2).We choose Then the system (2) takes the form (after some simplification) where The rest of the paper is structured as follows.Section 3 is composed of a blend of dynamical behaviours of the model (3) in deterministic environment.In Section 4, computer simulations are carried out to validate some of the analytical results.The stochastic version of the model is formulated in Section 5 to take into account the effect of fluctuating environment characterized by Gaussian white noises and the analysis of the stochastic model by the method of statistical linearization (Valsakumar et al. [51]).The criteria for non-equilibrium fluctuation and stability are also derived in this section.Section 6 contains the general discussions of the paper and a comparative study of the stability in the deterministic and stochastic environment.

Dynamical behaviour
In this section, we present a blend of dynamical behaviours of the model (3).Haque and Li [41] have derived some important results for the model (3) (although they have omitted some proofs for the lack of space).We have also derived some useful results.The results are listed below.

Theorem 1. All solutions of the system (3) which start in R 2
+ are uniformly bounded.Haque and Li [41] have omitted the proof of this theorem for the lack of space.We have provided the proof in the Appendix.
Theorem 2. The system (3) always has the boundary equilibrium points E 0 (0, 0) and E 1 (1, 0).The interior equilibrium point E * (X * , Y * ) exists uniquely if and only if the following condition is true: In this case X * and Y * are given by

Theorem 3. The necessary condition to reach the origin following
See the proof in Haque and Li [41].
The proof of the theorem is deferred to the Appendix.
Haque and Li [41] have also obtained the above result on permanence but they have omitted the proof.Readers can have the proof in the Appendix.Theorem 6.If E * exists, then it is locally asymptotically stable or unstable according as The proof is given in the Appendix.The proof is deferred to the Appendix.
Remark.When a Hopf bifurcation occurs, there exists small amplitude periodic orbits near E * (X * , Y * ).

Numerical simulation
In this section we present computer simulation of some solutions of the deterministic system (3).We choose the parameters of the system (3) as δ = 1.5, β = 1.8, γ = 0.2.Then E * (X * , Y * ) = (0.2800, 0.1867) and ∆ = 0.03000 > 0. Therefore, by Theorem 6, E * is locally asymptotically stable.The phase portrait for different choices of X(0) and Y (0) is depicted in Fig. 1.Clearly the solution is a stable spiral converging to E * .If we gradually increase the value of β (keeping other parameter values fixed), it is observed that the behaviour of the system changes as β passes through the bifurcation value β * = 1.875 (obtained by using Theorem 8).For β = 1.9 > β * , we see that ∆ = −0.0200< 0 and therefore by Theorem 3, E * (X * , Y * ) = (0.2480, 0.1653) is unstable.The corresponding phase portrait is a periodic orbit around E * (see Fig. 2).Thus using β as control, it is possible to break the stable behaviour of the system and drive it to an unstable state.Also it is possible to keep the populations at a desired level using the above control.

The stochastic model
It is now well understood that deterministic models in ecology do not usually incorporate environmental fluctuations based on the idea that in the case of large populations, stochastic deviations (or the effect of random environmental fluctuation) are small enough to be ignored.A stochastic model provides a more realistic picture of a natural system than its deterministic counterpart.We have already mentioned that environmental fluctuations play an important role in the interaction of M. occidentalis -T.mcdanieli (Wollkind and Logan [44], Wollkind et al. [42], Collings [43]).Svirezhev and Logofet [52], Dimentberg [33], Samanta [32], Maiti and Samanta [39], and many others have mentioned that the basic mechanism and factors of population growth like the resources and vital rates -birth, death, immigration and emigration, change non-deterministically due to random environment.Here we assume that fluctuations in the environment will manifest themselves mainly as fluctuations in the growth coefficients of the prey and predator since these are the main parameters subject to coupling of a prey-predator pair with its environment (Svirezhev and Logofet [52], Dimentberg [33]).Thus, as a starting point of this section, we consider the following stochastic version of the model (3): Clearly the system (4) has the same interior equilibrium point as the system (3).To study the behaviour of the system (4) about the steady state E * (X * , Y * ) we put X = u+X * and Y = v+Y * .Then the system of equations ( 4) reduces to the form (neglecting power greater than 2 and in terms of the deviated variables u, v) (Bandypoadhyay and Chakrabarti [38], Maiti and Samanta [39]): where η 1 (t), η 2 (t) are assumed to be independent Gaussian white noises satisfying the conditions: Here ǫ j (j = 1, 2) are the intensities or strengths of the random perturbations, δ, the Dirac delta function and .represents the ensemble average and These are Itô type of stochastic differential equations (non-linear coupled bivariate Langevin equations) governing the system behaviour about the steady state E * (X * , Y * ).The solutions (u(t), v(t)) of ( 5) subject to known initial values (u(0), v(0)) represent the state of the system at time t > 0. Now, we are concerned with stochastic differential equations ( 5) which are driven by Gaussian white-noises and interpreted mathematically as Itô stochastic differential equations.Gaussian white noise, which is a delta-correlated random process, is very irregular and as such it is to be treated with care.In spite of this it is an immensely useful concept to model rapidly fluctuating phenomenon.Of course, true white noise does not occur in nature.However, as can be seen by studying their spectra, thermal noise in electrical resistance, the force acting on a Brownian particle and climate fluctuations, disregarding the periodicities of astronomical origin etc. are white to a very good approximation.These examples support the usefulness of the white-noise idealization in applications to natural systems.Furthermore, it can be proved that the process (u, v), a solution of (5), is Markovian if and only if the external noises are white.These results explain the importance and appeal of the white noise idealization (Horsthemke and Lefever [53]).
In the past few decades, different techniques of linearization of nonlinear stochastic differential equations giving rise to a set of deterministic moment equations have been receiving a great deal of attention in different fields of science and technology (Nisbet and Gurney [31], Haken [54]).Jumarie [55] pointed out the fact that moment techniques can be used to solve a large class of problems in stochastic optimization involved with the problem of stochastic optimal control.In the following, the behaviour of the stochastic model ( 4) about the steady state will be cultured by the technique of statistical linearization developed by Valsakumar et al. [51].This approach has some limitations in their validity compared to the original non-linear stochastic differential equations have.However, this technique has some advantages in reducing the complexity of the solution of original non-linear equations without loss of information about the system too much.

Statistical linearization: moment equations
The statistical linearization of the equations ( 5) consists of replacing the equations by the system of linear equations: where the errors in the above linearization are given by The unknown co-efficients p i , q i , s i (i = 1, 2) of the equations ( 6) are determined from the minimization of the averages of the squares of the errors (7).We determine the unknown co-efficients by demanding that (Valsakumar et al. [51], Van Kampen [56], Bandyopadhyay and Chakrabarti [38]): Also we use the following expressions (Valsakumar et al. [51]): Then the expressions for p i , q i , s i (i = 1, 2) are given by The co-efficients are the functions of the parameters involved with the model system and also of the different moment involving u and v. Simple calculations lead to the system of equations of the first two moments: where we have used the relations Let us now assume that the system size expansion is valid such that the correlations ε 1 and ε 2 given by ( 9) decrease with the increase of the population size and they are assumed to be the order of the inverse of the population size N (Valsakumar et al. [51], Baishya and Chakrabarti [57], Bandyopadhyay and Chakrabarti [38]): Therefore, using the expressions ( 8), ( 9) and keeping the lowest order terms and replacing the averages u and v by their steady state values u = v = 0 (Nicolis and Prigogine [58]), we get the following reduced equations for second order moments: where D stands for the operator d dt .

Non-equilibrium fluctuation and stability analysis
Eliminating u 2 and v 2 from the equations of ( 11), we get the following third order linear ordinary differential equation in uv : The auxiliary equation of ( 12) is given by where Then the nature and structure of the roots of ( 13) will solely depend upon the quantities A and H (since 2A 3 − 3AB + C = 0).

Case 1. H < 0.
In this case the roots of ( 13) are given by The solutions of the system of linear equations (11) are then given by where L ij (i, j = 1, 2, 3), P 1 , P 2 are constants.When a 1 < 0, then A > 0 and consequently each of u 2 , v 2 , uv decreases to zero with increasing time (as c 2 < 0).So according to the criteria of stability in the sense of second order moments, E * (X * , Y * ) is stable when a 1 < 0. When either a 1 > 0 or A < 0 (or both), then E * is unstable as the second order moments diverge with increasing time.Now we notice that A > 0 ⇔ ∆ < 0. Therefore when a 1 < 0, the deterministic stability criterion (∆ < 0) of E * is satisfied and it is enough to guarantee the stability for the stochastic model (4).Also if E * is unstable for the deterministic system (3), then it is also so for the stochastic system (4).
In this case the roots of ( 13) are given by The solutions of the system of linear equations ( 11) are then given by where N ij (i, j = 1, 2, 3), Q 1 , Q 2 are constants.If a 1 < 0, then A > 0 and consequently the deterministic stability criterion (∆ < 0) is satisfied.In the stochastic environment, however, it is seen that each of the second order moments u 2 , v 2 , uv converges with increasing time whenever a 1 < 0 and A > √ 3H.Thus if a 1 < 0 and A > √ 3H, then the stochastic system is stable in the sense of second order moments .On the other hand, whenever A < √ 3H with a 1 < 0 then each of u 2 , v 2 , uv diverges with increasing time and hence the stochastic system becomes unstable, although it is stable in the deterministic environment.In all other cases also, the stochastic system is unstable.Some of our results in stochastic environment are illustrated through numerical simulation.When we choose δ = 1.5, β = 1.5, γ = 0.2, then H = −0.0875, a 1 = −0.0400< 0, A = 0.2400 > 0, and consequently E * is stable in the sense of second order moments (see Fig. 3).Fig. 4 shows unstable behaviour of E * when δ = 1.9, β = 1.5, γ = 0.2 (H = −0.0639< 0).On the other hand, if δ = 0.9, β = 0.9, γ = 0.02 then H = −0.0205, a 1 = −0.3019< 0, A = 0.32190 > 0, and E * is stable in stochastic environment (see Fig. 5).  5. Behaviour of the system (4) in the sense of second order moments when δ = 0.9, β = 0.9, γ = 0.02.

Discussion
In ecology, there are various concepts of stability and various biological phenomena are connected with them (Svirezhev and Logofet [52]).In the present paper, we are concerned with the deterministic and stochastic dynamical aspects of stability of a modified version of the Holling-Tanner model.The deterministic part consists of the results on the boundedness, persistence, permanence, stability and bifurcation of the system under positive initial population distribution.It is seen that if the growth rate of the prey is high but the growth rate of the predator surpasses it, then the persistence of the system is guaranteed.Also, using β as control, it is possible to break the stable (spiral) behaviour of the system (3) and drive it to an unstable (cyclic) state.Also it is possible to keep the levels of the populations at a required state using the above control.It is interesting to notice that the condition for local stability of E * (X * , Y * ) is independent of the carrying capacity K for the prey.So, it can be said that a change in the carrying capacity will not change the stability of E * (X * , Y * ).Hence, "paradox of enrichment" can not happen to this system.A very interesting observation on the deterministic extinction can be made from the result of Theorem 3. It indicates that if the consumption rate is high, then this will drive the prey population to extinction, and consequently, the predator will die (extinct) in starvation.Nowadays, almost all the developing countries are increasingly realizing the potential of the method of bio-control for exotic pests because of the long list of side effects of the chemical pesticides.In bio-control, both the areas of co-existence and coextinction are very important.
For the stochastic version of the model system, that is, for the model ( 4) under random perturbation, it is observed that the deterministic criteria of stability is no longer enough to guarantee the stability of the positive interior equilibrium E * (X * , Y * ).The deterministic and stochastic systems behave alike with respect to stability of E * when H < 0 and a 1 < 0. When H > 0 with a 1 < 0, then the deterministic criterion for stability A > 0 (or ∆ < 0) is satisfied but this criterion is not enough to determine the stability in the stochastic environment.In this case, the stability of the stochastic system (4) requires an additional condition A > √ 3H besides the deterministic stability condition.But if A < √ 3H, then the system (4) becomes unstable.Thus the stability of the system under random perturbation changes as A passes through the value √ 3H.Thus, to sum up, we have two comparative studies here: (i) the modified Holling-Tanner model versus the traditional one, and (ii) the stochastic model (4) against its deterministic counterpart.In the first case, it is seen that "paradox of enrichment" cannot happen to the modified system, whereas it can happen in the classical one.Also, if the prey-catching capacity is higher then the intrinsic growth rate of the predator only, then the traditional Holling-Tanner model leads to a total extinction, but for the modified model, it depends on the growth rates of both prey and predator.Analyzing the stability and bifurcation results of the Holling-Tanner model and the modified one, we may roughly say that the modified model is more stable than the traditional one (Haque and Li [41]).In the second case, when the stability results on the stochastic model ( 4) is compared with those of its deterministic counterpart (i.e., the model (3)), it is observed that the deterministically stable system remains stable under stochastic perturbation if certain conditions (viz., H < 0, a 1 < 0 or H > 0, a 1 < 0, A > √ 3H) are fulfilled.On the other hand, if H > 0, a 1 < 0, A < √ 3H, then the random perturbation has a destabilizing effect on the system.Thus, roughly speaking, stability and instability are consequences of stochastic perturbation of the model system under consideration.Such a conclusion is in good agreement with Prajneshu [59], Baishya and Chakrabarti [57], Samanta [32,35], Samanta and Maiti [37], Maiti and Samanta [26] and many others.

Appendix
Proof of Theorem 1.Let (X(t), Y (t)) be any solution of the system with positive initial conditions. Since by a standard comparison theorem, we have Then, Thus, all the solutions of the system (3) enter into the region Hence it is the region of attraction in this case, proving the theorem.
Hence, by definition, the theorem follows.
Proof of Theorem 6.The variational matrix at E * (X * , Y * ) is It is easy to see that the trace of V (E * ) is and its determinant det V (E * ) = γX * > 0.
The characteristic equation of V (E * ) is where P = −tr V (E * ) and Q = det V (E * ).
Since Q = det V (E * ) > 0, it is clear that E * is locally asymptotically stable or unstable according as P > or < 0.
Hence the theorem.
Hence all the conditions of the Hopf bifurcation theorem (Murray [46]) are satisfied and the theorem follows.