Analysis of a Predator-Prey Model: A Deterministic and Stochastic Approach

This paper investigates the deterministic and stochastic fluctuations of a predator-prey model. The predator is experienced in hunting two different prey simultaneously. Each prey has logistic growth in the absence of the predator. The rate of experience of the predator in hunting each prey is varied using a simulated dataset. The deterministic and stochastic nature of the dynamics of the system are investigated. Stability analysis is performed, using slight perturbation around the non-zero, interior equilibrium point, to determine where the system loses stability. The variation of the predatory experience parameter causes the system to experience Hopf bifurcations. These stability changes and the addition of stochastic noise are explored using time series graphs. The co-existence and extinction of the populations are affected over time . Citation: Addison LM (2017) Analysis of a Predator-Prey Model: A Deterministic and Stochastic Approach. J Biom Biostat 8: 359. doi: 10.4172/21556180.1000359


Introduction
The rise and decline of world populations and ecological species alike have increased the need for sustainable development. One way to achieve this is through the use of mathematical models which analyse these interactions. The original Lotka-Volterra Lotka [1], Volterra [2] model has been used for many years to provide an insight into the dynamics of predator-prey interactions over time. This model has since been modified to include logistic growth rates in prey populations, functional responses of predators to their prey, harvesting rates of prey in different habitats, and time delays in finding or consuming prey. These models are inherently deterministic in nature but within recent times, the notion of random fluctuations or noise has been incorporated to reflect real -life situations.
Deterministic models in ecology do not usually incorporate environmental fluctuation based upon the idea that in the case of large populations, stochastic deviations (or effects of random environmental fluctuation) are small enough to be ignored [3]. Researchers have mainly been interested in the dynamical consequences of population interactions, often ignoring environmental variability altogether [4]. However, populations are affected by environmental noise. These effects may be more noticeable when the population size is small. Noisy or stochastic effects can have large implications on the qualitative dynamics of the system. Generally stochastic effects may enhance, diminish or even completely change the dynamical behaviour of the system [5]. White noise is one such instance of noise, defined as the generalized meansquare derivative of the Wiener process, which is also known as Brownian motion. Many papers have explored this idea in different aspects. Oksendal [6] discussed the idea of allowing for randomness in some of the coefficients of a differential equation for a more realistic mathematical model. Mandal and Banerjee considered additive and/or multiplicative environmental noise terms which allow the deterministic model to be extended to a stochastic one.
It follows that stochastic differential equation (SDE) models play a prominent role in a range of application areas, including biology, chemistry, epidemiology, mechanics, microelectronics, economics, and finance [7]. According to Bandyopadhyay and Chattopadhyay [3], there are two main methods to develop the stochastic model corresponding to the existing deterministic one to study the effect of fluctuating environment: replacement of the environmental parameters by some random ones or, the addition of the random fluctuation directly to the prey and predator growth equations without changing parameters. This paper employs the second approach, which is also used by Baishya and Chakrabarti [8] and Bandyopadhyay and Chakrabarti [9]. Consequently, our model is a predator-prey system in which the predator has the ability to attack two species of prey simultaneously. This is affected by the rate of returns to the experience of the predator in hunting the specific prey. This rate is based on the concept of diminishing returns, first developed in 1767 by French Economist Turgot [10]. Although first applied to agriculture and the environment, the idea is popular economic law underlying the use of resources both ecologically and economically. The notion is that if other variables remain constant, the returns to profits can dimish, if production is increased beyond a certain point. Applied to this model, a more general approach is taken in the simulated dataset, where this rate may be increasing or decreasing between the prey species.
In addition, the rate of returns to experience can be compared to a similar parameter discussed in the predator-prey investment model in Brander and De Bettignies [11] where investors 'prey' on venture capital opportunities based on their experience in the industry over time. Few papers have observed this rate in a prey -predator model in a biological sense. The main objective of our paper is to pave the way for a bioeconomic approach to the model by exploring the changes in the determinsitic and stochastic dynamical behaviours of the predatorprey model for different values of this rate. These are explored in terms of co-existence, extinction, the stability of the equilibrium, existence of bifurcation and the addition of stochasticity. Therefore, this paper demonstrates the use of the prey-predator PREY: P 1 and P 2 represent the number of prey of two different species; PREDATOR: X represents the number of predators of one species, where P 1 , P 2 , X ≥ 0.
Parameters used in the model are all non-negative: 2): Interaction term between each prey species and predator; β i (i=1,2): Returns to experience of predator in catching this species of prey; ρ i (i=1,2): Natural growth rate of prey; For simplicity, we assume the prey species do not interact.

Existence of positive interior equilibrium point
In order to find the steady state solutions of the equations in the system, set each equation to zero. This gives: Let the interior equilibrium point be * * * 1 2 ( , , ) P P X . For simplicity, the superscripts are have been dropped in calculations. The analysis by Dubey and Upadhyay [12] in used here.
In order to obtain a positive value of X, the following inequalities must hold: Using equation (7), when P 2 → 0, thenP 1 →P 1f , that is Now P 1f is positive and real if Also, using equation (7), This equation does not produce a closed form solution for P 1g analytically, so it will be solved numerically. Also, using equation (8),  The equations (7) and (8) have a unique point of intersection (P 1 ,P 2 ) if Using known values of P 1 and P 2 , the value of X can be calculated using: We may thus write the following Lemma, resuming the use of the superscript,*, to depict the positive steady state values:  (7) and (8) and satisfy the inequalities in equations (5), (6), (10) and (16).

Stability analysis of interior equilibrium
The stability of the interior equilibrium point is discussed by examining the equilibrium point * * * where u,v and w are small perturbations about the equilibrium point. Assuming Taylor's Theorem, all terms are expanded about the equilibrium point, while neglecting higher order terms of u,v and w.

Hopf bifurcation analysis
Suppose that j 11 <0, j 22 <0 and j 33 <0. Then a 1 >0, a 2 >0 and the characteristic polynomial equation has two purely imaginary roots if and only if a 1 a 2 =a 3 for some value of β 1 say This equation has the three roots: These roots can be represented, in general as follows: Applying Hopf's Bifurcation Theorem from Marsden and Mckracken [13] to the original system in equation (1), it is important to verify the traversality condition: Substitution of λ 1 (β 1 )=p(β 1 )+iq(β 1 ) into equation (26), calculation of respective derivatives and separation of real and imaginary parts gives: Therefore, and λ 3 (β 1 )=-a 1 (β 1 ) ≠ 0. Hence, Hopf Bifurcations are obtained for the parameter β 1 and a similar result can be obtained for β 2 .We will confirm this using numerical simulations.

Numerical simulation results
Numerical simulations were performed using the MATCONT (version matcont3p4) package in MATLAB [14] to verify analytical results. Numerical analysis is used in conjunction with Hopf Bifurcation analysis to determine, using a range of parameter values, if the two prey and predator can coexist in terms a steady state solution (a node or a stable focus) or a stable oscillatory solution (limit cycle) [15]. The values of the parameters are chosen purely for simulation purposes: .8, ρ 2 =0.5, α 1 =0.001, α 2 =0.08, δ=0.1, K The parameters β 1 and β 2 representing the experience of the predator with respect to hunting prey 1 and prey 2 respectively are used as bifurcation parameters. The parameters β 1 and β 2, respectively are varied between 0.1 and 1.0 in increments of 0.1, and the values where Hopf bifurcations occur for β 1 and β 2 are recorded in Tables 1 and 2 respectively. At these values, the dynamics of the system change from stable to unstable and there is a possibility for extinction in one or more populations.
In Table 1, Hopf bifurcations for the system where β 2 is the bifurcation parameter only occur when β 1 is varied from 0.6 to 1.0, respectively. The same occurs in Table 2 when β 1 is the bifurcation parameter. Therefore, the system, in both instances, does not experience stable equilibrium for these values, but mimics an oscillatory solution. The system is unstable for these values. Figures 1 and 2 show an unstable and stable case for each Table respectively. The time series graphs show the co-existence of the three populations of two prey species and one predator as time progresses for each parameter set.
In Figure 1, the system is unstable for β 1 =1.85, β 2 =0.7, that is, returns to experience for predator with respect to hunting prey 1 is greater than than for prey 2. The population of predators is larger than both prey species initially, with prey 2 species being the least in size. After some time, (time, t, greater than 500), the predator population increases drastically, while the two prey species also increase, but never out perform the size of the predator species in this time span. While all three populations co-exist, it can be noted that the prey 2 species is smaller than prey 1 and at certain times it decreases close to extinction, and then increases.
A similar argument can be made with respect to Figure 2. Here, β 1 =0.7, β 2 =1.4 (returns to experience for predator with respect to hunting prey 2 is greater than than for prey 1). In this case, the prey 1 species is greater than both predator and prey 2 populations initially for this parameter set. This trend continues for all populations with no major changes in the initial values over time for this parameter set. It must be noted here that time has no specific units and may represent minutes, hours, days and so on depending on the biological system being under consideration.

The Stochastic Model
Stochastic perturbations are introduced into parameters in the original model in (1). There are different ways of constructing a stochastic differential equation model corresponding to an existing deterministic one, in order to study the effect of environmental fluctuations [4]. Stochastic perturbations of this type were successfully applied to different mathematical models by Bandyopadhyay and Chakrabarti [9], Beretta et al. [16], Carletti [17] and Mukhopadhyay and Bhattacharyya [18].
Typically, the Ito stochastic differential equation has the form: where the solution (Y t ; t>0) is an Ito process, r is the drift coefficient, s is the diffusion coefficient and W(t)=(W 1 (t),W 2 (t),…,W d (t)) is a d-dimensional process with independent Wiener processes. The latter term models noise in the environment, sometimes called Gaussian white noise (generalized derivative of Brownian motion) and is useful in representing random fluctuations. Application of the stochastic differential equations system to the original deterministic model produces the following result: Where σ i , i=1,2,3 are real constants representing the size of noise in the system due to the environment and dW i (t), i=1,2,3 represent independent, standard Wiener or Brownian motion processes.   The asymptotic stability of the positive equilibrium point E * of the stochastic system (31) is investigated and compared to the dynamics of the deterministic model in (1).

Stochastic numerical method
There are various methods to find the stability of an equilibrium point for a stochastic differential equation. In this model, this takes the form of white noise type stochastic perturbations of the variables P 1 ,P 2 and X proportionally distanced from the positive equilibrium point ( ) , , E P P X = [19]. The analytical method uses the idea that the stochastic model has the same equilibrium points of the original deterministic model [20]. The interest in this work is to observe the behavior of solutions for the stochastic system around the deterministic equilibrium values by taking small perturbations and varying the noise parameter, σ i .
It is difficult to calculate analytical solutions to the non-linear stochastic model. Hence, a discretization scheme may be used to numerically simulate trajectories for the stochastic differential equations. The simplest method of obtaining an approximate stochastic solution is the Euler-Maruyama (EM) Method [21]. This method approximates a sample path for each Stochastic Differential Equation (SDE) over an interval [0,T] which is divided into n sub-intervals of size Δt, where T t n ∆ = . The point (P 1 (t i ), P 2 (t i ), X(t i ))=(P 1i , P 1i , X i ) for t i =0,t,2t,…,T, where i=0,1,2,…,n. Full details of the numerical method are provided [22].    Table 2. The parameters reflecting the returns to experience for the predator in hunting and catching prey 1 and prey 2 species areβ 1 =0.7, β 2 =1.4. All populations co-exist with no extinction.
The EM method has strong order of convergence 0.5 and it converges to the Ito solution since it is the simplest strong Taylor approximation [23]. The method approximates each differential equation at time, t i using: Applying the EM method to the stochastic system in equation (31) for an initial point (P 10, P 20, X 0 )=(P 1 (0) , P 2 (0) , X 0 ) gives the following: (1 ) , where i=0,1,2,…,n, and η i has the distribution of a standard Normal random variable, that is, N(0,1) and existence, uniqueness and convergence of the numerical approximations are assumed to be satisfied. For simplicity, assume σ 1 =σ 2 =σ 3 =σ, that is, the value of noise is the same for each equation.

Numerical simulation results
The numerical simulations in the stochastic model use the parameters from the deterministic model which is close to the bifurcation point so the system starts to experience fluctuations: Δt=0.008 is used in the simulations, which are repeated 10000 times up to a time, t=80.
When the strength of environmental noise is virtually non-existent, that is, very close to zero, the system behaves like the deterministic model. It is evident that when the strength of the noise parameter is increased, the fluctuations of the sample paths also increase in an erratic manner as seen for the values of σ=0.05, 0.05, 0.5 and 1.4 shown in Figures 3-5, respectively. Figure 5 shows that the number of predators and prey experience extinction or are exhausted after an oscillation of large amplitude at time, t=40 for this particular simulation. This varies with other simulation runs but populations are exhausted at some time in the future with this level of noise. Figures 6 and 7 demonstrate the variations of noise when the prey populations are plotted against the predator population, respectively. These plots also re-iterate the idea that increasing environmental noise in the system cause the oscillations of population numbers to become unstable. In Figure 6, as time increases, the trajectories oscillate around the steady state value for prey 1 in an extremely erratic fashion in comparison to Figure 7 in which the pattern, with noise, is similar to that of the original oscillations for the deterministic model. However, when noise is increased to a value of 0.5, which is relatively high, both phase portraits experience completely erratic fluctuations, more so for prey 2 than prey 1.

Discussion
In this paper, the deterministic and stochastic features of a predatorprey model were examined. The analytical solution to the deterministic co-existence equilibrium point was found and stability analysis was performed on this point. Numerical analyses of Hopf Bifurcation points were performed with respect to parameters representing returns to predator experience since this is an important measure of the growth of the predator population. Theoretically, the greater the experience a predator has in hunting and catching a prey species, the faster the predator population should increase. This is similar to the economic idea that the larger the number of investment opportunities (prey) available to a Venture Capitalist (predator) in an industry, the greater his experience in investing in that industry [11]. However, this notion may not always be the case since noise factors such as drastic environmental changes, disease in prey or predator or even prey defense mechanisms can affect the deterministic approach. Hence the idea of positive or negative returns to predator experience in hunting prey, as well as stability variations in the system.
The numerical results were presented for a certain set of parameters and the time series graphs demonstrated the stability and oscillatory co-existence of prey and predator. There exists variations in the stability dynamics of the deterministic model with the variation of the bifurcation parameter, 'experience' of the predator in hunting the prey species. These stability ranges can assist biologists and ecologists alike in sustainability policies with respect to conservation of species or resources. The simulation seeks to showcase the importance of keeping parameters within certain ranges for co-existence of predator and prey populations in the deterministic case.
Dynamical systems are affected by the intensity of noise in the system. Continuous perturbations give rise the white noise, which has the distribution of a Gaussian random variable. In the stochastic model, the qualitative behavior of the system near the steady state values is of practical importance. Low levels of noise have the effect of a deterministic model, however, as noise is increased, the system oscillates with this random variation. As the value increases, these erratic oscillations begin to die down. The fluctuations allow for the examination of the co-existence of the venture capitalists and their Time series for two prey and one predator when the value of noise is 0.05

Conclusion
In this paper, the dynamical behavior of a three dimensional deterministic prey-predator model has been investigated using Stability and Hopf bifurcation analysis and numerical simulations. The deterministic model shows that the predator can co-exist with two prey species for which he has different levels of experience hunting. Otherwise, the unstable cases can alert ecologists and policy makers about the parameters which would cause the system to oscillate with the possibility of extinction for one or more species.
The variation of the parameter representing the returns to experience of the predator in hunting its prey provides valuable insight into the changes in stability of the system. The limitations in the availability , size of population, over-hunting defense mechanisms of prey, as well as noise, can be affect the returns to the predator's experience in hunting prey. Predators may continue to hunt at the same rate but obtain less prey due to the aforementioned factors. Hence, this supports the study of noise in the model. Also, the introduction of noise can cause the system to experience sharp increases, decreases or extinction over short spaces of time. It is useful to study the qualitative behavior of the systems near the interior equilibrium point in order to allow for optimal strategies for managing the co-existence of predators with the available prey species. It also allows for the proper management and monitoring of parameter values, which can cause extinction. Value of noise = 0.5 Figure 6: Simulation of the effect of varying noise from low to high, that is, σ=0 to σ=0.5, on the phase portraits for predator plotted against prey 1 populations. Low noise shows oscillation along a straight line for the prey 1 species with respect to predator and as noise increases, the oscillations fluctuate more. Value of noise = 0.5 Figure 7: Simulation of the effect of varying noise from low to high, that is, σ=0 to σ=0.5, on the phase portraits for predator plotted against prey 2 populations. Low noise shows oscillation along a limit cycle for the prey 2 species with respect to predator, and as noise increases, the oscillations fluctuate more.