Stochastic analysis of a prey – predator model with herd behaviour of prey

In nature, a number of populations live in groups. As a result when predators attack such a population, the interaction occur only at the outer surface of the herd. Again, every model in biology, being concerned with a subsystem of the real world, should include the effect of random fluctuating environment. In this paper, we study a prey–predator model in deterministic and stochastic environment. The social activity of the prey population has been incorporated by using the square root of prey density in the functional response. A brief analysis of the deterministic model, including the stability of equilibrium points, is presented. In random environment, the birth rate of prey species and death rate of predator species are perturbed by Gaussian white noises. We have used the method of statistical linearization to study the stability and non-equilibrium fluctuation of the populations in stochastic model. Numerical computations carried out to illustrate the analytical findings. The biological implications of analytical and numerical findings are discussed critically.


Introduction
About an entire century has already been elapsed on understanding and analyzing the basic rule between living food and its eater.It is not possible to construct a mathematical model that will fit entirely any natural subsystem.To find out the most suitable model that may analyze and forecast the natural phenomena, ecological and mathematical researchers are performing continuous numerous study.This process was initiated by Alfred James Lotka [20] and Vito Volterra [32].The Lotka-Volterra prey-predator model has been consequently used as a machine to introduce numerous theoretical and mathematical concepts in population modelling.This model has been modified incorporating Verhulst or Logistic growth to take into account the fact that the resource is limited, and also to avoid the structural unstable nature of the model.
The mass action predation term (in Lotka-Volterra model), though reasonable, is not the only possibility.Depending upon the behaviour of populations, more suitable 'functional response' has been developed as a quantification of the relative responsiveness of the predation rate to change in prey density at various populations of prey.In this connection, Holling family of functional responses [13,14,15] are the most focused, and in particular, the Holling type-II functional response has served as basis for a very large literature in prey-predator theory [1,3,4,22,26,30].According to Holling, the probability of a given predator encountering prey in a fixed time interval T t depends linearly on the prey density.We may express this relationship in the form X = αT s X, where X is the number of prey consumed by one predator, X is the prey density, T s is the time available for searching and α is the constant of proportion, termed as 'search efficiency'.Assuming that, each predator requires a handling time T h for each individual prey that is consumed, the time available for searching T s is given by: T s = T t − T h X.And then the above expression leads to X = αT t X − αT h X X, which implies X = αXT t /(1 + T h αX).Thus, the functional response would be of the following form: where X(T ) is the prey density at time T , α is the search efficiency of the predator for prey, T h is the handling time for each prey.Handling time is the average time spent by the predator on processing the food, beginning from the time the predator finds the prey item to the time the prey item is eaten.The type-II functional response also includes the fact that a single individual can feed only until the stomach is not full, and so a saturation function would be better to describe the intake of food.
In natural ecosystems, many living beings live forming herds and all members of a group do not interact at a time.There are many reasons for this herd behaviour, such as searching for food resources, defending the predators, etc.As a consequence, it is necessary to search for suitable form of functional response to describe this social behaviour.Only a few works have so far tried to enlighten this area.These works have demonstrated an ingenious idea that suitable powers of the state variables can account for the social behaviour of the populations.For example, to explore the consequence of forming spatial group of fixed shape by predators, Cosner et al. [9] have introduced the idea that the square root of the predator variable is to be used in the function describing the encounter rate in two-dimensional systems.Similarly, for three-dimensional systems, the two-third power of the predator in the encounter rate would better describe such group behaviour by predators.Unfortunately, such an idea has not been used by the researchers for about a decade.The work of Chattopadhyay et al. [8] may be regarded as a strong recognition of this concept.Then came the most innovative works of Ajraldi et al. [2] and Braza [7], which gave such modelling a new dimension.We recall their central ideas in the next paragraph.
http://www.mii.lt/NALet X be the density of a population that gathers in herds, and suppose that herd occupies an area A. The number of individuals staying at outermost positions in the herd is proportional to the length of the perimeter of the patch, where the herd is located.Clearly, its length is proportional to

√
A. Since X is distributed over a two-dimensional domain,

√
X would therefore count the individuals at the edge of the patch.Thus, when attack of a predator on this population is to be modelled, the functional response should be in terms of square root of prey population.This is the main idea of Ajraldi et al. [2].Braza [7] has placed a strong emphasis on this concept, and he has introduced a new functional response, where the prey density in (1) is replaced by its square root.That is, the functional response takes the form In this paper,we have considered the original model suggested by Braza [7] as our basic deterministic model.Most of the mathematical models proposed and analyzed in ecological literature are in the unvarying framework assuming that the observed dynamics are driven exclusively by internal deterministic mechanism.But real biological systems will always be exposed to influences that are not completely deterministic.Hence, there is an increasing need to expand the deterministic system to models including stochastic influence or noise that embrace more complex variations in the dynamics.Mathematical modelling in randomly fluctuating environment elicited modellers to introduce the stochastic differential equations.The deterministic ordinary differential equations are thus extended to stochastic differential equations.This is done in two ways: (i) when relevant parameters are taken as suitable stochastic processes, or (ii) when the deterministic system are partly driven by noises.In this paper, we follow the second one.
In many cases, mainly in nonlinear dynamical systems the noise may drastically change the characteristic of the systems.Generally stochastic effects may enhance, diminish or even completely change the dynamical behaviour of the system.Lorenz [19] and May [24,25] are the initiators in this context.Renshaw [29] 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.
A central obstacle in stochastic modelling of ecosystems is that there is an insufficiency for proper mathematical tools to analyze the nonlinear multi-dimensional stochastic processes.However some different techniques of linearization of non-linear stochastic differential equations have been developed by some researchers (see [11,28,31]).Though these approaches have some limitations in their validity compared to the original nonlinear stochastic differential equation, these have been used to reduce the complexity of the solution of original non-linear stochastic differential equations without changing the characteristics of the system too much.
In this paper, we have considered a prey-predator model, where the 'functional response' is of the form Holling type-II, but the prey density is replaced by its square root.
The paper is organized as follows.In the next section, we have discussed the details of the assumptions in our basic deterministic model and the significance of the parameters used in it.The positivity and boundedness, stability of equilibrium points have also been carried out.Section 3 deals with the stochastic version of the model, which takes into account the effect of fluctuating environment characterized by Gussion white noises.The method for linearization developed by Valsakumar et al. [31] has been used to find out the criterion for stability.Section 4 contains the numerical verification of the analytical findings.Biological significance of the analytical and numerical findings have been discussed in Section 5.

The deterministic model
The model we consider here is originally proposed by Braza [7].However, in his paper, Braza has presented the analysis of a simpler form of the model, assuming the average handling time to be zero.In the following, we discuss how the model is constructed.
Let X(T ) denotes the prey population density at time T .In the absence of predators, the prey population is assumed to have a logistic growth with intrinsic growth rate r and environmental carrying capacity K. Let Y (T ) denotes the density of the predator that has the only food source X.Our main consideration is that the prey population live in herd.Therefore, we consider the model under the framework of the following pair of nonlinear differential equations: Here α is the search efficiency of the predator for prey, T h is the average handling time for each prey, c is the biomass conversion rate of prey population to predator population, and δ is the natural death rate for predator population.All parameters are assumed to be positive.
To reduce the number of parameters and to make it easier, we non-dimensionalize system (2) with the following scaling: Then system (2) takes the form (after some simplification) http://www.mii.lt/NAwhere

Positivity and boundedness
Theorem 1.All solutions of system (3) that start in R 2 + remain positive forever.The proof is simple, and therefore, it is omitted.The following theorem ensures the boundedness of system (3).
Theorem 2. All solutions of system (3) that start in R 2 + are uniformly bounded.Proof.Let (x(t), y(t)) be any solution of system (3).Since we have lim sup t→∞ x(t) 1. Hence, Applying a theorem on differential inequalities [6], we obtain and for t → ∞, 0 W γ d .
Thus, all the solutions of (3) enter into the region Hence, the theorem.

Equilibrium points and criteria for their existence
The following theorem states all the possible equilibrium points of system (3).
Theorem 4. System (3) always has the trivial equilibrium E 0 (0, 0) and axial equilibrium , the interior equilibrium point E * (x * , y * ) exists and is given by

The trivial equilibrium
At the origin E 0 (0, 0), the variational matrix becomes indeterminate.System (3) being not linearizable, due to the square root term, the stability of the origin cannot be evaluated.
Rescaling the variable x = p 2 , the singularity may be overcome [2].But such rescaling may hide the true dynamics in case of prey-predator system.
We think the study of Braza [7] is more realistic to highlight the effect of the square root term by a local nonlinear analysis of system (3) to uncover the singular dynamics near origin.As the population densities approaches origin, it is reasonable to assume x(t) sufficiently small with the initial value x 0 = x(0) near to origin so that (i) x 2 or higher order terms vanishes, (ii) 1 + a √ x ≈ 1 and (iii) b √ xy is negligible to y.Under this approach, system (3) shifted to If the prey population is considerably smaller than predator population, i.e. x = O(y α ) with α 2, then the prey population first extincts and the predator population follows suit.For α < 2, origin becomes a saddle causing system (3) to be unstable near origin.

The predator-free (i.e. prey-only) equilibrium
At E 1 (1, 0), the variational matrix V (E 1 ) is given by The corresponding eigenvalues are −1, −(d − b/(1 + a)), and hence, system (3) is stable at Thus, there is a threshold value of the predator mortality rate over which the predator population goes to extinction. http://www.mii.lt/NA

The interior equilibrium
At the interior equilibrium E * (x * , y * ), the variational matrix V (E * ) is given by where a 11 , a 12 and a 21 are given by Then we have the following theorem.
Proof.The characteristic equation of V (E * ) is given by where For existence of E * , we have d < b/(1+a) that ensures B > 0. And for the proposed condition, it is obvious that A > 0. Thus, the roots of (4), given by are negative or with negative real part.Hence, the theorem.

Stochastic model
It is impossible to think of the real world without environmental fluctuations, and hence, it is essential to investigate the influence of environmental noise in biological systems through proper modelling.It is well recognized that the basic mechanisms and factors, such as birth, death, etc., change non-deterministically due to environmental fluctuations.From this viewpoint, we have perturbed the prey growth rate and predator death rate by noises [10,21,22,23,33].Then the deterministic system (3) is modified to a stochastic system under the framework of the following stochastic differential equations: where the perturbed terms η 1 (t), η 2 (t) are assumed as the independent Gaussian white noises that satisfy the conditions η j (t) = 0 and η j (t)η j (t ) = j δ(t − t ) for j = 1, 2.
Here j (j = 1, 2) are the intensities or strength of the random perturbation, δ is the Dirac delta function and • represents the ensemble average.Now we are concerned with stochastic differential equations, 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 (x, y), 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 [16].
Using the transformation x = x * e u , y = y * e v (where (x * , y * ) is the equilibrium of the deterministic system (3)) and neglecting the terms with degree more than two, we have the following pair of stochastic differential equations in Itô form (non-linear coupled bivariate Langevin equation): where http://www.mii.lt/NA

Statistical linearization: moment equations
To study the complexity of non-linear stochastic differential equations, different techniques has been derived by several researchers and mathematicians [12,28].The mostly used way is to linearize it, which gives rise a set of deterministic moment equations.Jumarie [17] 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 ( 5) about the steady state will be cultured by the technique of statistical linearization developed by Valsakumar et al. [31].Though this is not the only way, but 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.Thus, we have the following linearized system: where the errors for linearization are given by The unknown coefficients p i , q i , s i (i = 1, 2) of equations ( 7) are determined from the minimization of the averages of the squares of errors (8).To compute the unknown coefficients we use the similar techniques used by Valsakumar et al. [31], Van Kampen [18], Bandyopadhyay and Chakrabarti [5]: Also using the following expressions for higher moments derived by Valsakumar et al. [31]: Nonlinear Anal.Model.Control, 21(3):345-361 the expressions for p i , q i , s i (i = 1, 2) are given by The coefficients are the functions of the parameters involved with the model system and also of the different moments involving u and v. Simple calculations lead to the system of equations of the first two moments where the following relations have been used: Let us now assume that the system size expansion is valid such that the correlations ε 1 and ε 2 given by ( 12) decrease with the increase of the population size and they are assumed to be the order of the inverse of the population size N , i.e.
Therefore, using the expressions ( 10), (11) and keeping the lowest order terms and replacing the averages u and v by their steady state values [27] u = v = 0, we get the following reduced equations for second order moments: http://www.mii.lt/NA

Non-equilibrium fluctuation and stability analysis
The characteristic equation of the coefficient matrix of system ( 13) is given by where To solve equation ( 14), let us replace µ by m − A and that leads to the following standard cubic:

Now two possible cases arise depending on the value of H.
Case 1: H < 0. In this case, all roots of equation ( 15) are real and given by 0, ± √ −3H.Hence, the roots of the original equation ( 14) becomes Here all the eigenvalues are real.They are negative if A > √ −3H > 0. Hence, system ( 13) is stable if A > 0 and A > √ −3H.Now A > 0 is the criteria for deterministic stability of E * , but in stochastic environment, it is not enough to ensure the stability of E * .So stability of the deterministic model does not guarantee the stability of the stochastic model.For example, if √ −3H > A > 0, the system is stable for deterministic environment but not in stochastic environment.On the other hand, if the system is unstable in deterministic environment (i.e.A < 0), the system must be unstable in random environment also.Case 1: H > 0. In this case, roots of equation ( 15) are given by 0, ±i √ 3H, and hence, the roots of the original equation ( 14) becomes Here the eigenvalues have negative real parts if and only if A > 0. Consequently, system (13) is stable in stochastic environment for A > 0. Further, since A > 0 is the condition for stability of the deterministic model; the system behaves alike (with respect to stability) in deterministic and stochastic environment.

Numerical simulation
Beside analytical findings, numerical simulations are also important; because simulations can be used to validate the analytical findings.For various choices of the parameters of the model, we have performed the simulations using MATLAB.It is observed that they are in good agreement with our analytical findings.
If we take the parameters as a = 0.2, b = 0.85, d = 0.5; then by Theorem 4, the interior equilibrium E * (x * = 0.4444, y * = 0.4198) exists as d is less than b/(a + 1) = 0.7083.It is also noticed that (ad + 3b)d 2 − (b + ad)(b − ad) 2 = 0.1281, a positive quantity; and so the condition of Theorem 5 is also satisfied.Therefore, for this choice of parameters, E * is locally asymptotically stable.The corresponding phase portrait is presented in Fig. 1a, which is clearly a stable spiral converging to E * .Figure 1b shows that x and y approach their equilibrium values in finite time.For the above choices of parameters and including a random perturbation with mean zero, the behaviours of x and y with time are depicted in Fig. 2.This figure shows that both the prey and predator populations trace random paths with mean x * and y * .The stability of the stochastic model for these parameters are also established as the second order moments approach zero (see Fig. 3).Here the computed value of H = 0.1574 > 0 and A = 0.1340 > 0 satisfy the stochastic stability criterion.
Changing the parameters a little bit as (a = 0.2, b = 0.95, d = 0.5 with x(0) = 0.4, y(0) = 0.5) computed value of the expression (ad + 3b)d 2 − (b + ad)(b − ad) 2 becomes negative (= −0.0211), and hence, E * becomes unstable and there is a periodic orbit near E * (see Fig. 4a).The oscillations of x and y with time is shown in Fig. 4b.For these parameters, as H = 0.1950 > 0, the stability of the stochastic system is similar to that of the deterministic system.The unstable nature of the stochastic system has been shown in

Concluding remarks
So called 'functional response' is the key component in prey-predator ecosystems.It determines the basic characteristics of a system since both the prey death rate and predator growth rate are maintained by this.Depending upon the nature of populations in concerned real subsystem, various sophisticated functional responses are proposed and discussed in ecological literature throughout the last century.But the social activity like herd behaviour of the prey population for defense or food searching purpose has not been studied considerably.Here we have studied the dynamics of a prey-predator model (where the prey has a herd behaviour) in both deterministic and stochastic environment.To describe the herd mechanism, we have taken a modified Holling type-II functional response, where the prey density is replaced by its square root.This functional response is derived from the concept that the prey-predator interactions take place only at the outer surface of the herd formed by the prey.
The deterministic part consists of the basic findings like boundedness and stability of the system under positive initial population distribution.It is shown that all solutions of the system are uniformly bounded, which, in turn, implies that the system is biologically well behaved.It has been shown that (Theorem 3) the predator population dies off if d > (b/a); i.e. in terms of original parameters, δ > (c/T h ).Hence, the predator population is washed out from the system if its death rate is greater than the ratio of the biomass conversion rate and the prey handling time.It is obvious that due to disease or any other natural circumstances, if predator death rate increases significantly; whole predator population may wipe out.But the predator population may also go to extinction if its biomass conversion rate becomes low due to digestive problem or the predator takes much time in handling the prey.Stability at the origin is an important feature of http://www.mii.lt/NA the 'square root functional response'.At the origin, the Jacobian being indeterminate, eigenvalues cannot be computed to discuss the stability at origin.Here we follow the pathway of approximation suggested by Braza [7].We assume that near the origin, both the populations being small, their products or higher order terms may be neglected.And this gives an ecologically sound result that, if the prey density becomes smaller compared to the predator (i.e.x = O(y α ) with α 2), the system approaches its trivial equilibrium through stable path.This is a basic difference from the other models, where origin being a saddle, no such possibility arises.Here the prey population first go to extinction and then the predator population follows suit (see [7]).Stability of the deterministic model at the predator-free and interior equilibrium has been discussed critically.It is again seen that if the death rate of the predator population becomes greater than some threshold value (d > b/(1 + a), which is automatically satisfied if d > (b/a)), the system stabilizes at the predator-free equilibrium.On the other hand, A > 0 ensures the deterministic stability of the interior equilibrium (see Fig. 1).
In fluctuating environment, we have formulated the stochastic version of the model by including the additive Gaussian white noises that perturb the growth rate of the prey population and the death rate of the predator population.Then the resulting model is cultured by the method of statistical linearization developed by Valsakumar et al. [31].This leads to a system of linear differential equation in terms of second order moments.The criteria for stochastic stability of E * (x * .y* ) is derived.It is observed that if H > 0, the populations in deterministic and stochastic environments behave alike with respect to stability.However, if H < 0, the deterministic stability criterion (A > 0) is not sufficient to guarantee the stability in the stochastic environment, rather, a further restriction A > √ −3H is necessary for stochastic stability.Thus, if √ −3H > A > 0, the deterministically stable system become unstable under stochastic perturbation.Hence, to keep ecological balance in a fluctuating environment, the populations of the system are to be regulated in such a way that A > 0 when H > 0 and A > √ −3H when H < 0. These analytical findings have been verified by computer simulation (see Figs. 1-6).

Figure 6 .
Figure 6.For a = 0.2, b = 0.95, d = 0.5, unstable behaviour of the stochastic system (5) in the sense of second order moments.

Fig. 5 .
Fig.5.In Fig.6, the second order moments oscillate with increasing amplitude indicating instability of the stochastic model.