# The Holling Type II Population Model Subjected to Rapid Random Attacks of Predator

Authors: Jevgeņijs Carkovs; Jolanta Goldšteine; Kārlis Šadurskis
Journal: Journal of Applied Mathematics
 (2018)

---

## Abstract

We present the analysis of a mathematical model of the dynamics of interacting predator and prey populations with the Holling type random trophic function under the assumption of random time interval passage between predator attacks on prey. We propose a stochastic approximation algorithm for quantitative analysis of the above model based on the probabilistic limit theorem. If the predators’ gains and the time intervals between predator attacks are sufficiently small, our proposed method allows us to derive an approximative average dynamical system for mathematical expectations of population dynamics and the stochastic Ito differential equation for the random deviations from the average motion. Assuming that the averaged dynamical system is the classic Holling type II population model with asymptotically stable limit cycle, we prove that the dynamics of stochastic model may be approximated with a two-dimensional Gaussian Markov process with unboundedly increasing variances.


---

## Body

## 1. Introduction

One of the most popular models of the dynamics of interacting predator and prey populations, such as those found within invertebrate and similar domains of life, in mathematical biology is the system of ordinary differential equations proposed in [1]:(1)dxdt=rx1-K-1x-mxA+xy,dydt=-δ+γmxA+xy,where phase variables x and y denote the density of prey and predator populations, respectively. In this model it is assumed that in the absence of a predator the prey population has a potential carrying capacity K and develops according to the logistic law with an intrinsic growth rate r, and in the absence of a prey the predator population exponentially decreases to zero with the intrinsic growth rate δ. The mutual influence of changes in the densities of the prey and predator in model (1) is considered by the trophic function mxy(A+x)-1, where the positive parameter m is the prey consumption rate by the predator or, in other words, corresponds to the number of prey individuals that can be “eaten” per unit of time. The positive parameter A reflects the saturation of the amount of prey consumed and, in addition, depends on the rate of reaction of the predator, i.e., the time between attacks on the prey. The parameter γ in formula (1) is a conversion factor that determines the effect of “eaten prey” on the growth rate of the population of the predator. The popularity of the model (1) is explained by the fact that under certain assumptions about positiveness of parameters r,K,m,A,δ, and γ it is structurally stable [2]; that is, there exists a unique asymptotically stable periodic trajectory. This model describes stable fluctuations of the size of the predator and prey populations that are sometimes observed in biological ecosystems. In accordance with the continuous type dynamic system (1), both populations are in permanent contact, and the benefits gained and losses suffered by predator and prey correspondingly in an arbitrary small time interval [t,t+Δ) are proportional to the length of the interval Δ. In fact, it is clear that changes in the size of both populations are accidental and can only be modelled on average by formula (1). Therefore, subsequently papers were published that took into account the random nature of the model under study by adding stochastic terms of the white noise type in the system of (1) (see the review in [3]). It is important to note that the choice of this type of random perturbation allows preserving the principle of predictability of the behaviour of populations, since stochastic differential equations define random Markov type processes. In the case of modelling by the means of stochastic differential equations, the predator’s gain and the prey’s loss during time Δ contain not only terms proportional to Δ but also terms that are proportional to the increment of the Brownian motion process Δw≔w(t+Δ)-w(t). In most of these papers, the authors manage to prove the possibility of the existence of a stable stationary close to periodic ergodic process that describes the behaviour of the stochastic model under sufficiently small random perturbations.In the modern literature on mathematical biology many authors use stochastic differential equations as a mathematical model for predator-prey ecological systems (we refer again to the review in [3]). The most typical papers using stochastic models of predator-prey populations are [4–16]. The authors of these papers study the effect of stochastic perturbations of various parameters of classical models, adding either white noise to a chosen parameter [4–10] or the integral over a centered Poisson measure [11–16]. Using the apparatus of modern stochastic analysis, authors study the possibility of the existence of positive solutions, the stability of possible stationary solutions, and the existence of moments of solutions, as well as estimating the asymptotics of the moments of solutions and other properties of solutions. It should be noted that most of the aforementioned papers deal with models with stochastic additives containing no higher than second-order phase variables. If, however, the diffusion coefficients or the integrand of the integral over the Poisson measure has a more complicated form, then the apparatus proposed in the aforementioned papers can hardly be used. At the same time, when analysing the dynamics of some predator-prey biological communities, it may be necessary to investigate the consequences of nonpermanent random contacts that occur at random time moments. In this case, it is natural to assume that at the time of the contact extraction of prey by predator will be a random process that is proportional to the functional response and depends on the phase coordinates in a more complex form.However, in this kind of model, the predator’s gain and the prey’s loss during timeΔ are proportional to the normally distributed random variable Δw(t) with parameter Δ and therefore can be either positive or negative, which is poorly consistent with the definition of these terms in a formula of type (1). In this paper we propose a model that also makes it possible to take into account the stochastic character of the trophic function of the dynamics of populations as considered in the Holling type II model, but the predator’s gain and the prey’s loss are limited positive random variables. Besides, our model takes into account the possibility that the predator may take some time to attack the prey, and therefore there are intervals when both populations develop independently. Moreover, our model also fulfils the Markov property.We propose a method of approximate numerical analysis of the probabilistic characteristics of population dynamics, based on application of the asymptotic diffusion approximation algorithm [17]. Application of the diffusion approximation method for the asymptotic estimation of the probability characteristics of the Markovian dynamical system, analogously to the algorithm of the classical central limit theorem, consists of two steps. Initially, using the small parameterepsilon and the limit theorem for sequences of Markov processes, we find a deterministic dynamical system for the averaged phase variables. This is followed by finding a stochastic dynamical system for normalized deviations from solutions of the averaged dynamical system. The resulting stochastic differential equation of the Ornstein-Uhlenbeck type is well studied and may be relatively simply analysed [17].


## 2. The Model

Let us propose a model where the contacts between predator and prey occur very often at random time moments{0=τ0<τ1<τ2<⋯<τk<⋯}, which are the moments of discontinuity of the trajectories of a stationary piecewise constant Poisson process {ξ(t),t≥0} [18] given on a certain probability space (Ω,F,P) with an exponentially distributed length of the constancy intervals P(Δk-1>t)=e-ε-1t, where Δk-1=τk-τk-1, and ε is a small positive parameter. Let us denote {Ft,t≥0} as the minimal filtration [19] to which the process {ξ(t),t≥0} is adapted. Let us also assume that at time moments {τk,k∈N} the random variables {ξ(τk),k∈N} have a continuous distribution P(ξ(τk)≤ξ)=F(ξ) and, therefore, without any loss of generality, we may consider that this distribution is uniform on the interval [0,1]. It is known [18] that the probability distributions of a homogeneous Markov process {ξ(t),t≥0} are uniquely determined by its weak infinitesimal operator given by(2)Qεvξ≔limt↓01tEvξt-vξ∣ξ0=ξ.In our case the Markov process is given by an infinitesimal operator(3)Qεvξ=ε-1∫01vz-vξdz,where ε is a positive parameter.Let us proceed to the formal definition of the Markov dynamical system dealt with in this paper. Let us suppose that the time of observation starts att0≥0, and densities of prey and predator populations at this moment are x(t0),y(t0), respectively. Suppose that the size of the prey’s population {xε(t),t∈[t0,t0+τ1)} changes in accordance with the logistic equation(4)dxεtdt=rxεt1-K-1xεt.And the size of the predator’s population {yε(t),t≥0} changes by the rule(5)dyεtdt=-δ+γhξtxεtA+xεtyεt,where ξ(t)~R(0,1), h(ξ(t))>0, E{h(ξ(t)}≡m, and E{h2(ξ(t))}≡σ2. At the time moment t1=t0+τ1 the collision between predator and prey individuals occurs, and predator’s gain is given by the expression (6)xεt1-xεt1-=-εyεt1hξt1-xεt1-A+xεt1-.Here and further f(t-)≔lims↑0⁡ft+s. Up to the time moment t2=t1+τ2, the dynamics of the size of the prey and predator populations are also given by (5)-(6) with initial conditions x(t1) and y(t1). At the moment t2=t1+τ2 the predator finds the prey again and its gain is given by expression (6) where the argument t1 is replaced by t2. Then up to the moment t3=t2+τ3 the sizes of predator and prey populations change again according to law (4)-(5), etc.It is useful to note that, using terminology and results of monograph [20], Markov dynamical system (3)-(4)-(5)-(6) may be defined by the system of stochastic differential equations(7)dxεt=rxεt1-K-1xεtdt-ε∫01yεthuxεtA+xεt-1νεdt,du,dyεt=-δ+γxεtA+xεt-1∫01huduyεtdt,where νε(dt,du) is a Poisson random measure with the parameter ε-1dudt.Figure1 shows the trajectory of the solution of the system of the impulse-differential equations (4)-(5)-(6) with the initial conditions xε0=20 and yε(0)=3 and the parameter values ε=0.01, Ehξτk=m=5,  A=18,  γ=2,  δ=4,  r=1, and K=1. Figure 2 shows the trajectory of the solutions of the system of (1) with initial conditions x(0)=20 and y(0)=3 and with the same the parameter values m=5,  A=18,  γ=2, and  δ=4.Figure 1
Phase trajectory of Markov process (4)-(5)-(6).Figure 2
Limit cycle of dynamical system (1).As shown in the web-publication [21] and as it can be seen in Figure 2, the system of (1) with the aforementioned parameter values has a unique asymptotically stable periodic trajectory. In Figure 1 the trajectory of the system of (4)-(5)-(6) corresponding to the same initial conditions over time is also grouped within some neighbourhood of the limit cycle described above. Now let us show that the trajectories of the system of (4)-(5)-(6) for positive sufficiently small values of the parameter ε are located within a close neighbourhood of the trajectories of system (1) and coincide in average with them. Let us estimate the deviations from the averaged trajectories. The two-dimensional stochastic process {xε(t),yε(t),t≥0} is filtration {Ft,t≥0} adapted, and for any t>0,s>0,a∈R and b∈R we have the identity:(8)Pxεt+s>a,yεt+s>b∣xεt-u,yεt-u,u≥0≡Pxεt+s>a,yεt+s>b∣xεt,yεt;i.e., {xε(t),yε(t),t≥0} has the Markov property [18]. In addition, the dynamic characteristics (4)-(5)-(6) of the process {xε(t),yε(t),t≥0} are homogeneous in time. Consequently, the dynamical system defined by the Poisson process {ξ(t)} and (4)-(5)-(6) defines a homogeneous two-dimensional Markov process with the infinitesimal operator [18] given by(9)Lvx,y≔limt↓0⁡1tEvxεt,yεt∣xε0=x,yε0=y-vx,y.In a sufficiently small time t interval the stationary Poisson process ξ(t) can perform jumps from point ξ(0)=η to point ξ1, and both random variables are independent and equally uniformly distributed on the segment [0,1]. Therefore, when t↓0 the following asymptotic equalities hold:(10)Pξt≤x∣ξ0=η,t<Δ1=Pη≤xPt<Δ1=Pη≤xe-ε-1t=Pη≤x+ot,Pξt≤x∣ξ0=η,t≥Δ1=Pξ1≤x1-e-ε-1t+ot=ε-1tPξ1≤x+ot,where o(t) is an infinitely small quantity with order higher t when t→0. Consequently, if v(x,y) is a finite sufficiently smooth function, then(11)Evxεt,yεt∣xε0=x,yε0=y-vx,y=Evxεt,yεt∣xε0=x,yε0=y,Δ1≤tPΔ1≤t-vx,y+Evxεt,yεt∣xε0=x,yε0=y,Δ1>tPΔ1>t-vx,y=Evx+trx1-K-1x-εyhξ1xA+x+ot,y+t-δ+γhξ1xA+xy+ot-vx,y1-e-tε-1+ot+Evx+trx1-K-1x+ot,y+t-δ+γhξ1xA+xy+ot-vx,ye-tε-1+ot;substituting these asymptotic equations into formula (9), one can obtain an expression for a weak infinitesimal operator for a fixed ε>0:(12)Lεvx,y≔limt↓0⁡1tEvxεt,yεt∣xε0=x,yε0=y-vx,y=ε-1Evx-εyhξ1xA+x-1,y-vx,y+vx′x,yrx1-K-1x+vy′x,y-δ+γxEhξ1A+x-1y=ε-1∫01vx-εyhξxA+x-1,ydξ-vx,y+vx′x,yrx1-K-1x+vy′x,y-δ+γmxA+x-1y.


## 3. Asymptotic Analysis of Population Dynamics

As previously mentioned, we assume that the values of the time intervals{Δk-1,k∈N}, which the predator spends on search for prey, have an exponential distribution with a parameter ε-1; that is, E{Δk-1}=ε. If we assume that these intervals are relatively small, that is, the parameter ε is sufficiently small, then it is possible to apply the asymptotic methods proposed in [22] for analysing the Markov process {xε(t),yε(t),t≥0} given by the infinitesimal operator (12). At first it is necessary to check whether such a limit operator L0 exists, that, for any continuously differentiable function v(x,y) in some vicinity U0 of each point from the positive quadrant x≥0,y≥0, the following equality is true:(13)limε→0supx,y∈U0Lεvx,y-L0vx,y=0.In our case, this limit exists and it is equal to(14)L0vx,y=limε→0Lεvx,y=vx′x,yrx1-K-1x+vy′x,yy-δ+γmxA+x+vx′x,yymxA+x.Operator (14) corresponds to a deterministic dynamical system of differential equations(15)dx-dt=rx-1-K-1x--mx-A+x-y-,dy-dt=-δ+γmx-A+x-y-coinciding with the system of (1). Consequently [22], it can be asserted that for any positive number T there exist such positive numbers ε0 and M that for all ε<ε0 the following inequality exists:(16)supt<Tε-1Exεt-x-t+yεt-y-t≤Mε.This means that for averaged analysis of the behaviour of the studied population over a sufficiently large time interval one can use the system of (1), i.e., the results given in [2]. As shown in this paper, provided that A(mγ+δ)<K(mγ-δ), the phase portrait of the system of (15) is from the stable limit cycle, the orbit of which all other trajectories approach asymptotically. If the initial data {xε(0),yε(0)} for the system of (4)-(5)-(6) are chosen sufficiently close to the coordinates {x-,y-} of the aforementioned cycle, then the line {xε(t),yε(t),t≤Tε-1}, which corresponds to the solutions of these equations, does not leave a sufficiently small vicinity of this cycle. It was shown in [22] that the deviation of the real values {xε(t),yε(t)} from the mean values {x-(t),y-(t)} for each {t∈(0,Tε-1)} has the order of smallness ε and can be approximated by the solution of the corresponding system of linear stochastic differential equations. To construct these equations, new phase variables are introduced:(17)Xεt=xεt-x-tε,Yεt=yεt-y-tε.By definition (17) a two-dimensional random process {Xε(t),Yε(t)} on each interval t∈(τk-1,τk),k∈N, is given by the equations(18)dXεdt=1εrxε1-K-1xε-rx-1-K-1x-+mx-A+x-y-,dYεdt=1ε-δyε+γhξτk-1xεA+xεyε+δy--γmx-A+x-y-but at the moments of contacts τk,k∈N, between a predator and a prey it is given by the condition of a jump(19)Xετk-Xετk-=1εxετk-xετk-=-1εyετkhξτkxετk-A+xετk-.Substituting in formulas (18)-(19) the expression {xε(t),yε(t)} with the expression(20)xεt=Xεtε+x-t,yεt=Yεtε+y-twe can obtain the following system of equations:(21)ddtXεt=1εrXεtε-2rK-1x-tXεtε-εrK-1Xε2t+mx-tA+x-ty-tdYεtdt=1ε-δYεtε+y-t+γhξτk-1Xεtε+x-A+Xεtε+x-tYεtε+y-t+δy-t-γmx-tA+x-ty-tXετk-Xετk-=-εYετkε+y-τkhξτkεXετk-+x-τkA+εXετk-+x-τk.The process {x-(t),y-(t),Xε(t),Yε(t)} depends only on its state at the time moment t; that is, by definition (17) the two-dimensional random process {Xε(t),Yε(t)} is filtration {Ft,t≥0} adapted and, for a given solution of the system of (15), has the Markov property. Henceforward, it is convenient to study the behaviour of solutions of the system of (21) together with the system of (15), investigating the dynamic properties of the multidimensional Markov process {x-(t),y-(t),Xε(t),Yε(t)}. A weak infinitesimal operator L(ε) [18] of this process can be found by calculating the limit expression(22)LεVx,y,X,Y≔limt↓01tEVx-t,y-t,Xεt,Yεtx-0=x,y-0=y,Xε0=X,Yε0=Y-vx,y,X,Yfor an arbitrary finite sufficiently smooth function V(x,y,X,Y). Considering the possibility of contact between populations in a sufficiently short time interval (0,t), the following expressions can be used to calculate the limit (22):(23)EVx-t,y-t,Xεt,Yεtx-0=x,y-0=y,Xε0=X,Yε0=Y-Vx,y,X,Y=EVx-t,y-t,Xεt,Yεtx-0=x,y-0=y,Xε0=X,Yε0=Y∣Δ0≤t-Vx,y,X,Y+EVx-t,y-t,Xεt,Yεtx-0=x,y-0=y,Xε0=X,Yε0=Y∣Δ0>t-Vx,y,X,Y=EVx+trx1-K-1x-mxA+xy+ot,y+t-δ+γmxA+xy+ot,X+t1εrXε-2rK-1xXε-εrK-1X2+mxA+xy-εYε+yhξ1εX+xA+εX+x+ot,Y+t1ε-δYε+y+γhξ1Xε+xA+Xε+xYε+y+δy-γmxA+xy+ot-Vx,y,X,Y1-e-tε-1+ot+EVx+trx1-K-1x-mxA+xy+ot,y+t-δ+γmxA+xy+ot,X+t1εrXε-2rK-1xXε-εrK-1X2+mxA+xy+ot,Y+t1ε-δYε+y+γhξ1Xε+xA+Xε+xYε+y+δy-γmxA+xy+ot-Vx,y,X,Ye-tε-1+ot+εQVx,y,X,Y.Besides,(24)EεX+xaΔ0+εX+x-1=ε-1∫0∞εX+xat+εX+x-1e-tε-1dt=∫0∞εX+xaεs+εX+x-1e-sds=xA+x-1+Oε.Substituting the asymptotic by t expansions (23)-(24) into expression (22) after passing to the limit, one can obtain the formula for the weak infinitesimal operator of the homogeneous Markov process {x-(t),y-(t),Xε(t),Yε(t)}:(25)LεVx,y,X,Y,ξ=rx1-K-1x-mxA+xy∂∂xVx,y,X,Y+-δ+γmxA+xy∂∂yVx,y,X,Y+1εEVx,y,X-εhξ1Yε+yεX+xA+εX+x,Y-Vx,y,X,Y+1εrXε-2rK-1xXε-εrK-1X2+mxA+xy∂∂XVx,y,X,Y+1ε-δYε-γmxA+xy+γmYε+yεX+xA+εX+x∂∂YVx,y,X,Y.Since we are discussing the behaviour of populations for sufficiently small ε, then it is necessary to use a passage to the limit in the previous formula for tending the parameter ε to zero on the right. Therefore, one can use the smoothness of the function V(x,y,X,Y) and further use the asymptotic equalities:(26)EVx,y,X-εYε+yhξ1εX+xA+εX+x-1,Y=Vx,y,X,Y-εyxA+x-1m∂∂XVx,y,X,Y+ε-YxA+x-1m-yXA+x-1m+yxmXA+x-2∂∂XVx,y,X,Y+ε12yxA+x-1Eh2ξ1∂2∂X2Vx,y,X,Y=Vx,y,X,Y-εyxA+x-1m∂∂XVx,y,X,Y-εYx+AyXA+x-1mA+x-1∂∂XVx,y,X,Y+ε12yxA+x-1σ2∂2∂X2Vx,y,X,Y+oε.Using these expansions, formula (25) can be rewritten in the form convenient for passage to the limit:(27)LεVx,y,X,Y,ξ=rx1-K-1x-mxA+xy∂∂xVx,y,X,Y+-δ+γmxA+xy∂∂yVx,y,X,Y+rX-2rK-1xX-Yx+AyXA+x-1mA+x-1∂∂XVx,y,X,Y+12yxA+x-1σ2∂2∂X2Vx,y,X,Y+-δY+γYmxA+x+XγymAA+x2∂∂YVx,y,X,Y+Oε.It follows that the limit operator L=limε→0⁡L(ε) corresponding to the Markov process {x-(t),y-(t),X(t),Y(t),t≥0}, given by a system of ordinary differential equations (15), is a stochastic linear inhomogeneous stochastic differential equation(28)dXt=g11tXtdt+g12tYtdt+ftdwtand an ordinary homogeneous differential equation(29)dYtdt=g21tXt+g22tYt,where(30)g21t=γy-tmAA+x-t2,g22t=-δ+γmx-tA+x-t,g11t=1-2Kx-tr-Ay-tmA+x-t2,g12t=-x-tmA+x-t,ft=σy-tx-t2A+x-t,and w(t) is a standard Brownian process. Since the limit cycle is an asymptotically stable state of the system of equations (15), it is natural to assume that the points x-(t),y-(t),t≥0 in the system of equations (28)-(29) are the coordinates of the limit cycle mentioned previously. Using the results of [22], one can assert that on an interval with order 1/ε the probabilistic characteristics of the solutions xε(t),yε(t) of the impulse-differential equations system (4)-(5)-(6) with accuracy of the order of smallness greater than ε can be approximated by random processes x-(t)+εX(t),y-(t)+εY(t). In Figure 3 a sample trajectory of the process x-(t)+εX(t),y-(t)+εY(t),t∈[0,20] is simulated in MATLAB environment with the initial conditions x-(0)=20,  y-(0)=3,  X(0)=0,  Y(0)=0,  varX(0)=0,   varY(0)=0, covarX0,Y(0)=0,  and ε=0.01 on the phase plane {x,y}.Figure 3
Phase trajectory of Markov processx-(t)+εX(t),y-(t)+εY(t),t∈[0,20].Comparison between Figures1 and 3 indicates that the proposed approximation sufficiently reflects the dynamics of real processes. However, over a large time interval this approximation can differ significantly from the real process. This fact will be explained using the moments of the solutions of the system of (28)-(29). The probabilistic characteristics of the original process and the proposed approximation should be close [22] over a time interval with order 1/ε. This means that at each time moment t∈0,T/ε the centered second moments of the process xε(t),yε(t) can be approximated using the second moments of the process X(t),Y(t) in the form of asymptotic equalities:(31)varxεt=εq11t+oε3/2,varyεt=εq22t+oε3/2,covarxεt,yεt=εq12t+oε3/2q11t,where q11(t)=var{X(t)},q22(t)=var{Y(t)}, and q12(t)=covar{X(t),Y(t)}. For further analysis, using the notations(32)Gt=g11tg12tg21tg22t,Z→t=XtYt,φ→t=ft0,let us rewrite the system of (28)-(29) in a vector-matrix form:(33)dZ→t=GtZ→tdt+φ→tdwt.Further attention will be focused on the solution of the system of (28)-(29) with zero initial data. It is known [17] that this solution can be represented in the form of a stochastic Ito integral(34)Z→t=∫0tUtU-1sφ→sdws,where U(t) is the matrix solution of the ordinary differential equation(35)ddtUt=GtUtwith respect to the initial value U(0)=1001. By the definition of a stochastic integral [18], formula (34) defines a Gaussian process with zero mathematical expectation and a covariance matrix(36)Z→t=∫0tUtU-1sφ→sdws.The covariance matrix Q(t)=EZ→(t)Z→T(t)≔q11(t)q12(t)q12(t)q22(t) that corresponds to the integral (36) as a function from parameter t can be represented in the form of integral(37)Qt=∫0tUtU-1sφ→sφ→TsU-1sTUTtdsor in the form of a matrix solution of a system of ordinary differential equations(38)ddtQt=GtQt+QtGTt+φ→tφ→ttwith the initial condition Q(0)=0. Figures 4 and 5 show the time evolution of the variances q11(t)=E{X2(t)} and q22(t)=E{Y2(t)}.Figure 4
Trajectory of the varianceq11(t)=E{X2(t)}.Figure 5
Trajectory of the varianceq22(t)=E{Y2(t)}.Figure6 shows the time evolution of the covariance q12(t)=E{X(t)Y(t)}.Figure 6
Trajectory of the covarianceq12(t)=E{X(t)Y(t)}.Figures4–6 show that the variances of the normalized deviations increase very rapidly and even at t≥20 it cannot be recommended to use the most popular applied statistics “3σ rule”; that is, the inequality xεt-x-t≤3εq11t,yεt-y-t≤3εq22(t) because the quantities either x-(t)-3εq11(t) or y-(t)-3εq22(t) can become negative, which contradicts the definition of random processes xε(t) and yε(t). For clarity, Figure 7 shows the process q11(t),q22(t) trajectory on the time interval [0,30] on the phase plane {x,y}.Figure 7
Phase-plane trajectory for variances{q11(t),q22(t)}.Here the variable3εq11(t)=0.3q11(t) at some time point already exceeds the value 0.310000=30, although the value x-(t) (the abscissa of the limit cycle in Figure 2) does not exceed 24.Unbounded increase of variancesq11(t),q22(t) is a consequence of a resonance in the equation (38). Let ω be a period of the limit cycle x-(t),y-(t). Then U(t)UT(t) and φ→(t)φ→T(t) are 2ω-periodic matrix functions and their elements may be decomposed into a Fourier series by cosnπ/ωt and sinnπ/ωt,  n∈N. Therefore, the integrand in (37) contains a constant component and limt→∞⁡Q(t)=∞. At a finite time interval that depends on ε we can use the diffusion approximation (39)x-t+εXt,y-t+εYt,t∈0,Tεfor estimation of finite dimensional distributions of the process xε(t),yε(t).Even if a deterministic limit cycle closely approaches one of the axes (Figure8) we can use diffusion approximation x-(t)+εX(t),y-(t)+εY(t) (Figure 9) for prognoses of the population growth up to time T=20,ε=0.001.Figure 8
Phase portrait for (1),  m=1,r=1,A=0.5,K=1,γ=1, and δ=0.5.Figure 9
Phase trajectoryx-(t)+εX(t),y-(t)+εY(t),t∈[0,20], m=1,r=1,A=0.5,K=1,γ=1, and δ=0.2. x-(0)=0.57, and y-(0)=0.5.


## 4. Conclusions

During the analysis of population dynamics models, one must find a compromise between accuracy and complexity of the model. Numerical methods are helpful but their use actualises the stability issues of dynamical systems.In the case of interacting invertebrate type populations of predator and prey, the fact that the time between predator attacks on prey is random leads to a stochastic model. Even though the deterministic model has an asymptotically stable limit cycle, introducing random parameters into the model changes its trajectories significantly. The algorithm proposed and demonstrated in this paper allows one to analyse the average motion of a system as well as the random deviations from it and to find all characteristics of the resulting two-dimensional Gaussian Markov process. This process has an unboundedly increasing variance, leading to possible significant deviations even on relatively short time intervals.



---

*Source: 6146027-2018-06-20.xml*