Deterministic Chaos Versus Stochastic Oscillation in a Prey-Predator-Top Predator Model

The main objective of the present paper is to consider the dynamical analysis of a three dimensional prey-predator model within deterministic environment and the influence of environmental driving forces on the dynamics of the model system. For the deterministic model we have obtained the local asymptotic stability criteria of various equilibrium points and derived the condition for the existence of small amplitude periodic solution bifurcating from interior equilibrium point through Hopf bifurcation. We have obtained the parametric domain within which the model system exhibit chaotic oscillation and determined the route to chaos. Finally, we have shown that chaotic oscillation disappears in presence of environmental driving forces which actually affect the deterministic growth rates. These driving forces are unable to drive the system from a regime of deterministic chaos towards a stochastically stable situation. The stochastic stability results are discussed in terms of the stability of first and second order moments. Exhaustive numerical simulations are carried out to validate the analytical findings.


Introduction
"Ecological systems are open systems in which the interaction between the component parts is nonlinear and the interaction with the environment is noisy. This intrinsic nonlinearity can give rise to a complex behaviour of the system, which becomes very sensitive to initial conditions, various deterministic external perturbations, and to fluctuations present in nature" (Spagnolo et al., [36]). Mathematical modelling and nonlinear dynamics provide powerful tools in the analysis of such problems and many different models have been used to study the dynamics of interacting species in a variety of different contexts (see, e.g., May, [26]; Murray, [27]; Renshaw, [33]; Turchin, [40]). Many ecological interactions such as those between predator and prey, hosts and parasites, and hosts and parasitoids, are able to generate fluctuations. This is also the case in model systems whose rate of disturbance depends on local biological conditions, such as population densities (Possingham et al., [31]). On the other hand, in population dynamics of natural ecosystems, no clear example of deterministic chaos has been ever observed. Ellner and Turchin, [14] have developed a methodology and tools to detect deterministic chaos in a short and noisy time series. These authors analysed laboratory and field data extracted from published records and found that deterministic chaos is not a dominant explanation for the variability exhibited by these records, although many of these could support the phenomenon (Rai, [32]).
Predator-prey interactions often lead to oscillatory dynamics in temporal systems and in spatial systems when the rate of movement is large, so that individuals are effectively well-mixed and space becomes unimportant (Pascual et al., [29]). When individuals are not well mixed, however, properties of fluctuations in population densities, and in particular their amplitudes are known to vary with the spatial scale at which the model system is observed (Pascual et al., [29]). Modeling population dynamics in random environment is a way of studying the fluctuations of population size that has been affected by the stochasticity of external factors like weather (Abbas et al., [1]). Chaotically fluctuating populations face the danger of becoming extinct from two potential sources: demographic and environmental stochasticity. However, recent theoretical and experimental evidence suggest that chaos cannot be responsible for population extinction. Moreover, theoretical investigations showed that large random noise added by the environment seems to be a decisive factor responsible for the extinction of population (Lande, [21]; Allen et al., [3]). All natural systems are noisy. Chaotic dynamics in appearance is noise-like. The noisy character of these systems is caused either by intrinsic factor or extrinsic forces. When a system is subjected to external influences (deterministic or stochastic), complex dynamical phenomena (chaos and noisy chaos) emerge. These phenomena are often indistinguishable. This leads to failure of attempts to observe chaos (Upadhyay, [42]).
May, [26] analysed a predator-prey system under stochastic fluctuation considering white noise for population, and observed that a greater deviation from equilibrium shows instability (the deviations from mean population increases). The ecological effects for terrestrial and aquatic systems depend on the character of the physical frequency distributions as the general qualitative response of these systems could be inherently different. For terrestrial systems, the environmental variability is large at both short and long time periods and can be expected to develop internal mechanisms, which could cope with short-term variability and minimize the effects of long-term variations. Hence, analysis of the system with white noise gives better results. However for aquatic systems, less robust internal processes are needed to handle the smaller amplitude variability at short time periods commensurate with the life span of the organisms. Aquatic ecological systems are more likely to have coloured noise, instead of white noise as compared to terrestrial systems (Steel, [37]; Pimm, [30]). A population time series is said to be 'white' when no frequency dominates, 'red' when fluctuations are dominated by low-frequency, long term variations and 'blue' when high-frequency fluctuations dominate (Cohen,[13]; Lawton, [22]). There are two major mathematical schemes by which noise influences the dynamical systems, viz. output and dynamical noise. Such noise may take the form of an additive or multiplicative expression, illustrating the kind of parameters by which noise may appear in the equations of a dynamical system. In the case of an output noise, the geometrical structure of an attractor is stable which means that as the values of parameters increase, the dimension of the attractor also increases, but remains finite. On the other hand, in the case of dynamical noise, the geometrical structure of an attractor is unstable. This means that as the values of parameters are increased, the geometrical structure of the attractor vanishes (Argyris et al., [4]). Analyzes that explicitly taken into account the dynamical noise in the time series suggest that most ecological populations are not chaotic (Ellner and Turchin,[14]), but there have been very limited tests of these techniques with ecological models (Nychka et al., [28]). In this context, we will study the problem of existence of a stationary distribution of model system and its asymptotic stability.
In order to introduce our analysis we recall that earlier theoretical studies have shown that population dynamics in real ecosystems can be explained by the stochasticity arising from random fluctuations of environmental variables, and the complex and chaotic behaviour originating from the nonlinearities of the system. Many studies theoretically show that ecological factors such as habitat-heterogeneity, immigration and omnivore may impede and control chaos (Cai-lin and Zi-Zhen, [24]). Recent advances in modelling natural systems take into account for stochasticity in the dynamics of physical and biological variables, introducing noise sources which act directly on the species dynamics (multiplicative noise in the stochastic differential equations of the populations) and considering also random fluctuations in the dynamics of environmental parameters (Spagnolo et al., [36]). A number of investigations suggest that population dynamics is sensitive to noise colour (see, e.g., Xu and Li,[46]; Laakso et al., [20]). Based upon a few empirical records and simple forcing models Steele have suggested that terrestrial noise should be white, while marine noise should be red (Steele and Henderson, [38]).
Researchers have mainly been interested in the dynamical consequences of population interactions, often ignoring environmental variability altogether. However, the role of environmental fluctuations has recently been recognized in theoretical ecology. The deterministic models of various ecological systems do not incorporate environmental fluctuations, because in case of large populations, stochastic deviations are small enough to be ignored (Bandyopadhyay, [6]). Further, advances in community dynamics will require the use and analysis of nonlinear stochastic models, as illustrated by the studies reviewed by Bjornstad and Greenfell, [11]. The environmental fluctuations or noise, via their interaction with the nonlinearity of the system, has given rise to new counterintuitive phenomena: stochastic resonance (Valenti et al., [44]), noise-enhanced stability (Agudov and Spagnolo, [2]) noise-induced nonequilibrium transitions (Van Den Broeck et al., [45]), noise-induced multistability and nonequilibrium phase transitions (see, e.g., Li and Hanggi, [23]; Berry, [10]).
The main objective of the present paper is to investigate the existence of deterministic chaos and stochastic oscillations in a food chain consisting of three trophic levels and involving Holling type IV functional response. Organization of this paper is as follows: in Section 2 we describe briefly the deterministic model system and also show the boundedness of the model system. Local stability analysis and conditions for the Hopf bifurcation are derived in Section 3. In Section 4, we show the results of numerical simulations. In Section 5, we introduce the stochastic extension of our model. Significant outcomes and a comparative study between deterministic and stochastic model system are given in the concluding section.

Model System
We consider the following system as a model simulating a tritrophic food chain. The dynamics of a food chain consisting of three species and characterized by a Holling type-IV functional response, is governed by the following system of differential equations, where x(t) is the population density of the lowest trophic level species (prey) at time t, y(t) is the population density of the middle trophic level species (intermediate predator) at time t and z(t) is the population density of highest trophic level species (top predator) at time t. The intermediate predator y feeds on the prey x according to the Holling type-IV functional response and the top predator z preys on y according to the same functional response.
where a 1 , b 1 , a 2 , w, w 1 , w 2 , w 3 , c, i and a are positive constants. a 1 is the intrinsic growth rate of the prey population x, a 2 is the intrinsic death rate of predator population y in the absence of the only food x, the parameter c is the decay rate of the top predator z in the absence of its prey y and w 3 is the maximum value which per capita increase or growth rate of z can attain, and the ratio w 3 /w 2 is measure of its assimilation efficiency. w is the maximum value which per capita reduction rate of x can attain, w 1 has similar meaning to w with per capita increase rate. w 2 is the maximum value which per capita reduction rate of y can attain. b 1 is the strength of intra-specific competition among individuals of the prey species x. The parameter a can be interpreted as the half-saturation constant in the absence of any inhibitory effect. The parameter i, in turn, is a direct measure of the predator's immunity or tolerance of the prey. The model system (2.1)-(2.3) has ten parameters in total and it can be written in the following form: Obviously the right hand sides of the model system (2.4)-(2.6) are smooth functions on R 3 Therefore, the positive octant of the interior of R 3 + is an invariant region. Further the boundedness of system (2.4)-(2.6) is shown via the following theorem.
where p = min(a 2 , c). By Comparison Lemma we obtain, for all t ≥T ≥ 0, IfT = 0, then for all t > 0 we get Thus all the species are uniformly bounded for any initial value in R 3 + .
(ii) The equilibrium point E 1 = (a 1 /b 1 , 0, 0) exists on the boundary of the first octant.
(iii) E 2 = (x,ỹ, 0) is the planar equilibrium point on the xy-plane, wherex is the positive root of the equatioñ clearly E 2 exists provided the condition 0 <x < a 1 /b 1 is satisfied.
(iv) The nontrivial equilibrium E 3 = (x * , y * , z * ) exists if and only if there is a positive solution to the following set of equations Straight forward computation show that x * is a positive root of the cubic equation Now this equation can be written as there is a positive root of (3.6) that lies in (0, a 1 /b 1 ) when y * < aa 1 /w is satisfied. Also Therefore the positive equilibrium point E 3 = (x * , y * , z * ) exists uniquely under the following conditions Now in order to investigate the local behaviour of the model system (2.4)-(2.6) around each of the equilibrium points , the variational matrix V of the point (x, y, z) is computed as The eigenvalues of V 0 are a 1 , −a 2 , −c. There is an unstable manifold along the x-direction and a stable manifold along yz-direction. Therefore, the equilibrium point E 0 is a saddle point. The variational matrix for E 1 is given by From the variational matrix V 1 , it is found that the equilibrium point E 1 is locally asymptotically stable provided The variational matrix about the equilibrium point E 2 = (x,ỹ, 0) is given by The root of the characteristic equation p 3 (λ) = 0 of the above variational matrix about E 2 = (x,ỹ, 0) satisfy the following: The equilibrium point E 2 = (x,ỹ, 0) is stable or unstable in the positive direction orthogonal to the x − y plane, i.e. z -direction depending on whether λ 3 = w 3ỹ / ỹ 2 /i +ỹ + a − c is negative or positive, respectively. For the positive equilibrium point E 3 = (x * , y * , z * ) the variational matrix is given by According to the Routh-Hurwitz criterion, E 3 = (x * , y * , z * ) is locally asymptotically stable provided the following conditions are satisfied It is easy to show that, A 1 > 0, A 2 > 0, A 3 > 0 iff the following condition is satisfied By substituting the values of a ij , i, j = 1, 2, 3 we get (3.9) Clearly, M 1 > 0 provided that condition (3.7) is satisfied and M 2 > 0 provided that condition (3.8) is satisfied. Hence, if condition (3.7) and (3.8) hold, then the necessary and sufficient condition for (3.12) We summarize the above results in the following theorem, In order to investigate the Hopf bifurcation of the model system (2.4)-(2.6), we follow the technique given by Liu, [25]. According to this approach, the simple Hopf bifurcation at μ = μ * can occur provided A 1 (μ), A 3 (μ) and Ψ (μ) = A 1 (μ)A 2 (μ) − A 3 (μ) are smooth functions of μ in an open interval which includes μ * ∈ R such that Now, let the decay rate of the top predator c be the bifurcation parameter. If conditions (3.7) and (3.8) hold together with the condition We summarize the results for the existence of Hopf bifurcation in the following theorem:

Numerical Simulations
The model was numerically integrated to get the time series corresponding to the variables of model system (2.4)-(2.6) using Matlab 7.5. Since every nonlinear system has a finite amount of transients, the data points representing transient behaviour were discarded. Phase portraits were drawn using these data to obtain the geometry of the attractors. The chaotic attractor and its corresponding time series are depicted on the following data set (see Fig. 1 and Fig. 2): Via the bifurcation analysis of the model system (2.4)-(2.6), very rich and complex dynamics are observed, presenting various sequences of perioddoubling bifurcation leading to chaotic dynamics. For bifurcation diagram presented in Fig. 3, the value of the control parameter a 1 lies between 0.63 to 0.75, and the values of other parameter are the same as given in (4.1). Figure 3 shows clearly the evidence of the route to chaos through the cascade of perioddoubling. The period-doubling phenomenon leading to chaos is a well known    From the Hopf-bifurcation analysis results presented in Theorem 3, we observe that the dynamics of the model system (2.4)-(2.6) is also sensitive to the control parameter c, the death rate of the top predator z. Similar bifurcation diagram is drawn in Fig. 4. Here we note that the solution of the model system (2.4)-(2.6) is very sensitive to the death rate parameter c in the range 0.3 ≤ c ≤ 0.318. Moreover, it is evident that the route of chaos through the cascade of period doubling for decreasing value of the control parameter c and other important change in the chaotic set include an interior crisis in which a chaotic attractor undergoes a sudden increase in the size of attractor (Grebogi et al., [18]). The occurrence of such phenomena in iterated maps was discussed by Gottlieb, [17]. His algorithm for finding the cut-off parameter as boundary crisis helps us to conclude that the existence of chaos in the model system considered here is a result of crisis-limited chaotic dynamics. Closed curves marked as A and B in Fig. 4 correspond to invariant KAM tori (Upadhyay and Rai, [43]) in the phase space. When the bifurcation parameter c is decreased further, these curves break and give rise to chaotic dynamics (see, e.g., Upadhyay and Rai, [43]; Upadhyay, [41]).

Stochastic Analysis
There are various ways of constructing a stochastic differential equation model corresponding to an existing deterministic one, in order to study the effect of fluctuations in the environment. Researches studying environmental fluctuation in population dynamics adhere to two of these. The first way is to replace the environmental parameters in the deterministic model, by some random parameters (e.g. the growth rate parameter "r" can be replaced by r 0 + γ(ω, t), where r 0 is the average growth rate, γ(ω, t) is the noise function, and is the intensity of fluctuation). Secondly, one can add a randomly fluctuating driving force directly to the deterministic growth equations of prey and predator populations without altering any particular parameter (see, e.g., Baishya and Chakraborti [5]; Bandopadhyay and Chakraborti [7]; Tapaswi and Mukhopdhyay [39]). This kind of stochastic model formulation is known as noise-added/noise-driven system and this kind of modelling approach is capable of capturing the change in net growth rate of either population or change in different parametric values of the model, due to environmental driving force. A third kind of noise, multiplicative noise is also used sometimes, [19]. This is feasible when some of the coefficients in the problem are unknown apriori, and best modelled as random. The analysis here is often more involved than the additive case, which we will stick to in our setting.
We assume that stochastic perturbations of the state variables around their steady-state values E * are of Gaussian white noise type which is proportional to the distances of x, y from their steady-state values x * , y * respectively (Beretta et al., [9]) Similar stochastic perturbations are considered also in (Carletti,[12]). This type of stochastic perturbations firstly was proposed in (Beretta et al., [9]; Shaikhet, [35]) and was used later in (Bandyopadhyay and Chattopadhyay, [8]; Shaikhet, [34]; Carletti [12]) also in some other publications. The plus point of this method is that the equilibrium point x * , y * is a solution of stochastic system too.
We work on the following stochastic differential equation model system obtained from the deterministic model after introducing additive noise terms to the growth equations, to yield the following system of stochastic differential equation model subject to the non-negative initial condition x(0), y(0), z(0) ≥ 0 and σ i 's (i = 1, 2, 3), are intensities of environmental forcing. dw r (ω, t) = ξ r (ω, t) dt, r = 1, 2, 3, are path wise derivative of stationary Wiener processes. We suppress the ω dependence hence forth. Here ξ r (t) are mutually independent white noise terms characterized by ξ r (t) = 0 and ξ r (t)ξ s (t 1 ) = δ rs δ(t − t 1 ), (see, e.g., Gard [15]; Gardiner [16]). In order to minimize the mathematical expressions we rewrite the above system into the following matrix form and G = diag(σ 1 , σ 2 , σ 3 ). Here we have used x 1 , x 2 , x 3 for x, y, z, and f i are defined as follows In rest of this work we consider the stochastic differential system (5.4), whose solution (X t , t > 0) with non-negative initial condition X(0) = X 0 is a process on a probability space (Ω, F , P). F is a slowly varying continuous component or drift coefficient and G is the diffusion coefficient. Since the matrix G is independent of X, the system (5.4) is said to have additive noise (see, e.g., Bandyopadhyay and Chattopadhyay [5]; Carletti [12]).
In the previous section we have studied the sufficient condition for local asymptotic stability of coexisting equilibrium points and discussed the possibility of Hopf-bifurcating periodic solution, and finally the parametric domain within which system exhibits chaotic dynamics. Now we are in a position to determine whether the introduction of noise imposes greater irregularity in the dynamical behaviour, or if it contributes towards stabilization. It is worthy to note that the stochastic differential system (5.4) has no equilibrium point, and hence we consider a perturbation of the deterministic system (5.4) around the coexisting equilibrium point, and then add noise to these to formulate the governing differential equation for the first and second order moments for perturbation variables. Stability of trivial equilibrium for the system of differential equations will determine the amplitude of fluctuation of the populations around their deterministic steady state values in the presence of environmental fluctuation. Let (x * 1 , x * 2 , x * 3 ) denote the interior equilibrium point of the deterministic system whose components are given explicitly in earlier section. Substituting the perturbations 1, 2, 3), into the deterministic system, expanding it in Taylor series and omitting the second and higher order terms of small quantities we get the following system of equations: When noise is added to these perturbations we obtain, . One can easily verify that a 13 = 0 and a 31 = 0. Thus the above is derived via simple Taylor expansion followed by addition of noise, rather than a particular Ito or Stratanovich rule. From (5.8)-(5.10) we can write the system of ordinary differential equations for first order moments as follows Clearly the stability of interior equilibrium point for the deterministic system implies the stochastic stability of the model system (5.4) in terms of the first order moments. It is also interesting to note that the environmental driving forces have no role to determine the stochastic stability of the system if we introduce small perturbation from the equilibrium population density when the first order moments are concerned. This is clearly due to the fact that dw i = 0, i = 1, 2, 3 (5.11) Once we have the equations of the first order moments, in order to derive the governing differential equations for the second order moments, we use Ito's lemma to arrive at 14) du 1 u 2 = [a 21 u 2 1 + a 12 u 2 2 + (a 11 + a 22 )u 1 u 2 + a 23 u 1 u 3 ] dt + σ 1 u 2 dw 1 + σ 2 u 1 dw 2 , (5.15) du 1 u 3 =[a 32 u 1 u 2 +a 12 u 2 u 3 +(a 11 +a 33 )u 1 u 3 ] dt + σ 1 u 3 dw 1 + σ 3 u 1 dw 3 , (5.16) Since u i is a small perturbation, we make the additional assumption that This is physically realistic from the small magnitude of the perturbations, and will also help us in the forthcoming analysis. The above assumption is feasible via the following lemma.
Lemma 1. Consider the system of stochastic differential equations (5.12)-(5.17). Given a arbitarilly large but finite time T , there exists a constant C dependent only on T and the initial data, such that the following estimates hold uniformly for t ≤ T Proof. We demonstrate via Ito's lemma on the function F (u i ) = |u i | 2 2 , Integrating the above in the time interval [0, T ] yields Now adding up the above,taking expectations and using Cauchy Schwartz and Holder's inequality 1 yields Now using the integral version of Gronwall's lemma, and using the fact that Hence, via mean value theorem for integrals This completes the proof.
In order to proceed with the calculations of the second moment we demonstrate with (5.12). The analysis of the other equations is similar. We take expectations to yield d dt u 2 1 = 2a 11 u 2 1 + 2a 12 u 1 u 2 + σ 2 1 + 2σ 1 u 1 dw 1 .

(5.21)
Since u 1 and dw 1 are not independent, which is also true of the other variables involved, we cannot conclude that A more rigorous stochastic analysis is required if we want to show the above. We recap certain essentials to this end.  We now define a stopping time T for any given H as T = inf t : u 2 1 ≥ H . Next we integrate (5.21) from 0 to t ∧ T and take expectations to yield Of primary interest is the last term t∧T 0 Remark 1. M t in general is only a local martingale, and hence we cannot conclude M t = 0. In order to justify this we will first demonstrate that M t∧T is a uniformly integrable martingale, and then take a suitable limit via which M t∧T → M t . After which we will apply the earlier stated theory.
The first step follows via the following lemma.

Lemma 2.
The local martingale M t∧T , as defined via (5.24), is a uniformly integrable martingale.
Proof. We have that In order to derive the above given , we choose λ = C/ . Furthermore, the above follows via Chebyshev's inequality, Ito isometry and Lemma 1. This shows that M t∧T is a uniformly integrable martingale, and completes the proof.
We can now use the optional sampling theorem Theorem 4, to yield We now take H → ∞ and define a sequence of stopping times {T n } such that T n → ∞ and thus t ∧ T n → t. Incorporating these limits in (5.23) yield This implies that or equivalently, d dt u 2 1 = 2a 11 u 2 1 + 2a 12 u 1 u 2 + σ 2 1 . For the sake of the analysis we will assume u 1 (0) = 0. If this is not the case, then the initial condition can be incorporated into the σ 2 1 term. The analysis for the other equations proceeds as the above and we thus arrive at u 1 u 2 = a 21 u 2 1 + a 12 u 2 2 + (a 11 + a 22 ) u 1 u 2 + a 23 u 1 u 3 , d dt u 2 u 3 = a 32 u 2 2 + a 23 u 2 3 + (a 22 + a 33 ) u 2 u 3 + a 21 u 1 u 3 , d dt u 1 u 3 = a 32 u 1 u 2 + a 12 u 2 u 3 + (a 11 + a 33 ) u 1 u 3 .
Mathematically it is possible to find the solution of the above system of equations but it is quite impossible to write the analytical expressions componentwise. For notational convenience we assume that the steady-state for the second order moments is denoted by U * = (u 11 * , u 22 * , u 33 * , u 12 * , u 23 * , u 13 * ), where the steady state of u 2 r is denoted by u rr * (r = 1, 2, 3) and that of u r u s is denoted by u rs * , where r, s = 1, 2, 3 with r < s. Thus the stability of steady-states for second order moments depends solely upon the sign of real parts of eigenvalues of the matrix B defined by  Chaotic dynamics is observed in the ranges b 1 ∈ [0.68, 1.44]. Next we numerically integrate the stochastic model system with additive noise term for the same set of parameter values except b 1 = 2.3 using Euler-Maruyama scheme. We mention that co-existing equilibrium point E * (0.3275, 0.2133, 0.0327) is stable whenever b 1 = 2.3 and all other parameters are the same as mentioned above. Simulation results indicate that the distribution of population is not chaotic rather they fluctuate stochastically around some average value as shown in Fig. 5. Fig. 6 shows the distribution of prey, predator and top predator population over the range of values within which they fluctuate. It is clear that for stochastically driven system no steady state population distribution exists for either component or this happens for the reason that top-predator population distribution is low for most of the time compared to the prey and predators. Figure 6. Histogram for prey, predator and top predator population obtained from ten simulations run of the stochastic model system and data collected after neglecting the initial transient values.
Now we are in a position to interpret the analytical results with the help of the numerical simulations as well as by computing the matrix B and its eigenvalues. For the choice of parameters mentioned at the beginning of this paragraph except b 1 = 2.3, solving the system of equations (5.12)-(5.17) we find u 11 * = 0.00009345, u 22 * = 0.00004165, u 33 * = 0.00008638. Calculating the eigenvalues of B we find that two eigenvalues are positive and hence the second order moments are not stable and as a result the equilibrium population density obtained from deterministic analysis is not a stochastically stable equilibrium for the stochastic model system (5.1)-(5.3) in sense of the second order moments.

Discussions and Conclusions
In this paper, we have studied a three dimensional continuous time prey-predator-top predator model system, modelling a tritrophic food chain based on a Holling type-IV functional response. The stability analysis and Hopfbifurcation for non-negative equilibrium point has been analysed. The dynamical behaviour of the model system is investigated by using stability analysis, bifurcation theory and numerical simulation. Bifurcation analysis of model system (2.4)-(2.6) shows that, the system has rich dynamics including periodic and chaotic dynamics. Even if the selection of biologically realistic parameter values for the numerical simulation of ecological models is difficult and our parameter range is narrow, very rich and complex dynamics abound.
The effect of environmental fluctuation on the three species model system is observed after introducing additive noise terms to the growth equation of population. Numerical simulation of the stochastic model using Euler-Maruyama scheme suggest that the distribution of population is not chaotic rather they fluctuate stochastically around some average value. Thus stochasticity has important consequences and, depending on the intensity of noise, may overwhelm the stabilizing effects of control. If a chaotic population is to be controlled, it is necessary to know what kind of variability is involved, and consequently to estimate the strength and frequency of control necessary to override random effects.