Department of Mathematics and Statistics, Faculty of Science and Technology, The University of the West Indies, St. Augustine Campus, Trinidad, West Indies
Received Date: July 14, 2017; Accepted Date: July 25, 2017; Published Date: August 12, 2017
Citation: Addison LM (2017) Analysis of a Predator-Prey Model: A Deterministic and Stochastic Approach. J Biom Biostat 8: 359. doi: 10.4172/2155-6180.1000359
Copyright: © 2017 Addison LM. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Visit for more related articles at Journal of Biometrics & Biostatistics
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 .
Predator-prey; Stability; Deterministic; Stochastic; Noise
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 , Volterra  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 . Researchers have mainly been interested in the dynamical consequences of population interactions, often ignoring environmental variability altogether . 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 . 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  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 . According to Bandyopadhyay and Chattopadhyay , 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  and Bandyopadhyay and Chakrabarti .
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 . 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  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 predator – prey 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 model to investigate deterministic and stochastic dynamics. The stability of the system is presented analytically and Hopf bifurcations are investigated numerically where the system loses stability for particular values of the rate of experience of the predator. In addition, the dynamics of the system, with respect to the addition of white noise, are displayed and discussed. Simulations are performed and recommendations are made based on the results.
The following model is deterministic form of the model without noise:
PREY: P1 and P2 represent the number of prey of two different species;
PREDATOR: X represents the number of predators of one species, where P1, P2, X ≥ 0.
Parameters used in the model are all non-negative:
αi (i=1,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;
Ki (i=1,2): Maximum value of prey population i=1,2;
δ: Death rate of predator;
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 . For simplicity, the superscripts are have been dropped in calculations. The analysis by Dubey and Upadhyay  in used here.
In order to obtain a positive value of X, the following inequalities must hold:
Solving equations (2), (3) and (4) give
Using equation (7), when P2 → 0, thenP1→P1f, that is
Now P1f is positive and real if
Also, using equation (7),
Now if either:
1. A1>0 and B1>0, or
2. A1<0 and B1<0.
Similarly, using equation (8), when P2 → 0, then P1→P1g in the equation:
This equation does not produce a closed form solution for P1g analytically, so it will be solved numerically. Also, using equation (8),
Now if either:
1. A2<0 and B2<0, or
2. A2>0 and B2>0.
The equations (7) and (8) have a unique point of intersection (P1,P2) if
P1f <P1g. (16)
Using known values of P1 and P2, 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:
Lemma 4.1 The positive equilibrium point for real populations exists if and are both positive solutions of equations (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 , and are all positive. The equations are linearized using the substitutions:
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.
The characteristic equation has the form
For stability, it is necessary for the eigenvalues, λ in equation (19) to have negative real parts. These conditions are provided by the Routh- Hurwitz criteria, which states that a stable equilibrium occurs if and only if
a1>0, a3>0, a1a2 - a3> 0 (22)
We may thus write the following theorem:
Theorem 4.2 Given an equilibrium point of system (1), then once Lemma (2.1) holds and j11, j12, j13, j21, j22, j23, j31, j32, j33 are defined by equation (21), then the equilibrium point exists and is stable if and only if (22) holds where a1,a2 and a3 have been defined in (20).
Hopf bifurcation analysis
Suppose that j11<0, j22<0 and j33<0. Then a1>0, a2>0 and the characteristic polynomial equation has two purely imaginary roots if and only if a1a2=a3 for some value of β1 say . This value of β1 is unique. There exists an interval containing , say
for some ε>0 for which −ε > 0. This occurs at a2 > 0 for
. Hence, in a neighbourhood of , the characteristic polynomial cannot have positive real roots.
Now for , we have
This equation has the three roots:
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  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) . The values of the parameters are chosen purely for simulation purposes:
ρ1=0.8, ρ2=0.5, α1=0.001, α2=0.08, δ=0.1, K1=10, K2=10
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.
|β1||Stable(S) or Unstable(U) region/s||HB point/s for β1|
|0.1 to 0.6||Stable everywhere||No HB points|
|0.7||0<β1 ≤ 1.832349, (S)||1.832349,|
|1.832349<β1 ≤ 2.587247, (U)||2.587247,|
|β1 >2.704382, (S)|
|0.8||0 0<β1 ≤ 1.852958, (S)||1.852958,|
|1.852958<β1 ≤ 3.119409, (U)||3.119409|
|β1 >3.119409, (S)|
|0.9||0<β1 ≤ 1.875257, (S)||1.875257,|
|1.875257<β1 ≤ 3.475709, (U)||3.475709|
|1.0||0<β1 ≤ 1.892511, (S)||1.892511,|
|1.892511<β1 ≤ 3.82864, (U)||3.82864|
Table 1: Hopf Bifurcation (HB) point/s when β1 is the bifurcation parameter for values of different values of β2 respectively where ρ1=0.8, ρ2=0.5, α1=0.001, α2=0.08, δ=0.1, K1=10, K2=10.
|β1||Stable(S)/ Unstable(U) region||HB point/s for β2|
|0.1 to 0.5||Stable everywhere||No HB points|
|0.6||0<β2 ≤ 2.177352, (S)||2.177352,|
|2.177352<β2 ≤ 7.237154, (U)||7.237154|
|β2 >7.237154, (S)|
|0.7||0<β2 ≤ 2.179424, (S)||2.179424|
|2.179424<β2 ≤ 7.151218, (U)||7.151218|
|β2 >7.151218, (S)|
|0.8||0<β2 ≤ 2.178432,(S)||2.178432,|
|2.178432<β2 ≤ 7.066589, (U)||7.066589|
|β2 >7.066589, (S)|
|0.9||0<β2 ≤ 2.173640, (S)||2.173640,|
|2.173640<β2 ≤ 6.983508, (U)||6.983508|
|1||β2 >6.983508, (S)|
|0<β2 ≤ 2.164175, (S)||2.164175,|
|2.164175<β2 ≤ 6.902237, (U)||6.902237|
Table 2: Hopf Bifurcation (HB) points when β2 is the bifurcation parameter for values of different values of β1 respectively where ρ1=0.8, ρ2=0.5, α1=0.001, α2=0.08, δ=0.1, K1=10, K2=10.
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.
Figure 1: Time series for prey and predator populations for unstable case from Table 1. The parameters reflecting the returns to experience for the predator in hunting and catching prey 1 and prey 2 species areβ1=1.85, β2=0.7. All populations co-exist at certain points in time, while there are periods of extinction and growth for prey 2 species.
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.
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 . Stochastic perturbations of this type were successfully applied to different mathematical models by Bandyopadhyay and Chakrabarti , Beretta et al. , Carletti  and Mukhopadhyay and Bhattacharyya .
Typically, the Ito stochastic differential equation has the form:
where the solution (Yt; t>0) is an Ito process, r is the drift coefficient, s is the diffusion coefficient and W(t)=(W1 (t),W2(t),…,Wd(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 dWi(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 P1,P2 and X proportionally distanced from the positive equilibrium point . The analytical method uses the idea that the stochastic model has the same equilibrium points of the original deterministic model . 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 . 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 . The point (P1(ti), P2(ti), X(ti))=(P1i, P1i, Xi) for ti=0,t,2t,…,T, where i=0,1,2,…,n. Full details of the numerical method are provided .
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 . The method approximates each differential equation at time, ti using:
Applying the EM method to the stochastic system in equation (31) for an initial point (P10, P20, X0)=(P1(0), P2(0), X0) gives the following:
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:
ρ1=0.8, ρ2=0.5, δ=0.1, K1=K2=10, α1=0.001, α2=0.08, β1=0.7, β2=1.83.
The value of σ is varied between 0.05 and 1.4 in order to understand the effect of noise on the system of equations. The steady state values for the number of prey and predators fluctuate around the equilibrium point from the deterministic model. The time-step parameter.
Δ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.
Figure 3: Time series for the simulation of predator and prey populations where noise is low (σ=0.05). Populations fluctuate erratically with the addition of noise, but the system still progresses close to the original deterministic system. Prey 2 species experiences periods of extinction with time while prey 1 has the largest overall population size for this dataset. Predator species has steady with a sharp increase at t =50. All species continue to co-exist with time.
Figure 4: Time series for the simulation of predator and prey populations where noise is medium level (σ=0.5). The system begins to deviate more from the original deterministic system with the prey 2 species having sharp fluctuations in the population with time. The predator and prey 1 species both experience periods of extinction with time and periods of fluctuations. All species continue to co-exist with time.
Figure 5: Time series for the simulation of predator and prey populations where noise is high (σ=1.4). The original unstable case experiences few periods of drastic fluctuations for prey 1 species and prey 2. For this particular simulation, the predator and prey 2 experience early extinction over a long period of time. Prey 1 also experiences periods of drastic growth and decline and also extinction eventually. This level of noise promotes mass extintion in the populations.
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.
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.
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.
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 . 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 opportunities in the presence of noise. This noise can cause different periods of growth and decline of prey and predator species.
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.
Future work in this area can also explore the idea of different predatory functions for both the deterministic and stochastic models. The interaction parameters can also be converted to functions as well, to realistically mimic the evolution of an ecological species. The model can also be improved by incoporating interaction between the two prey species. Statistical techniques for parameter estimation can also be applied to these models and comparison of the model efficiency for different functional response rates is also possible.
By no means can the existing complex biological models be replaced. However, the use of non-linear differential equations and numerical analysis using deterministic and stochastic techniques with time series analysis can greatly assist in the study of the sustainability of species with a common predator.
Special mention and thanks to Professor Bhatt and Dr. Owen who have provided valuable feedback on this research. Also, special thanks to the Editors of the Journal of Biometrics & Biostatistics for this publication.