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

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.


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]: where phase variables  and  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  and develops according to the logistic law with an intrinsic growth rate , 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 ( + ) −1 , where the positive parameter  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  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 , , , , , 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 [,  + Δ) 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 Δ fl ( + Δ) − ().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][5][6][7][8][9][10][11][12][13][14][15][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][5][6][7][8][9][10] or the integral over a centered Poisson measure [11][12][13][14][15][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 Δ() 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 parameter epsilon 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].

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 < ⋅ ⋅ ⋅ <   < ⋅ ⋅ ⋅ }, which are the moments of discontinuity of the trajectories of a stationary piecewise constant Poisson process {(),  ≥ 0} [18] given on a certain probability space (Ω, F, P) with an exponentially distributed length of the constancy intervals (Δ −1 > ) =  − −1  , where Δ −1 =   −  −1 , and  is a small positive parameter.Let us denote {F  ,  ≥ 0} as the minimal filtration [19] to which the process {(),  ≥ 0} is adapted.Let us also assume that at time moments {  ,  ∈ N} the random variables {(  ),  ∈ N} have a continuous distribution ((  ) ≤ ) = () 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 {(),  ≥ 0} are uniquely determined by its weak infinitesimal operator given by In our case the Markov process is given by an infinitesimal operator 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 at  0 ≥ 0, and densities of prey and predator populations at this moment are ( 0 ), ( 0 ), respectively.Suppose that the size of the prey's population {  (),  ∈ [ 0 ,  0 + 1 )} changes in accordance with the logistic equation And the size of the predator's population {  (),  ≥ 0} changes by the rule Phase diagram with rapid random attacks, formulas ( 4)-( 6 4)-( 5)-( 6).
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 {  (),   (),  ≥ 0} is filtration {F  ,  ≥ 0} adapted, and for any  > 0,  > 0,  ∈  and  ∈  we have the identity: i.e., {  (),   (),  ≥ 0} has the Markov property [18].In addition, the dynamic characteristics (4)-( 5)-( 6) of the process {  (),   (),  ≥ 0} are homogeneous in time.Consequently, the dynamical system defined by the Poisson process {()} and ( 4)-( 5)-( 6) defines a homogeneous twodimensional Markov process with the infinitesimal operator [18] given by In a sufficiently small time  interval the stationary Poisson process () 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  ↓ 0 the following asymptotic equalities hold: = P ( ≤ ) +  () , where () is an infinitely small quantity with order higher  when  → 0. Consequently, if V(, ) is a finite sufficiently smooth function, then substituting these asymptotic equations into formula (9), one can obtain an expression for a weak infinitesimal operator for a fixed  > 0:

Asymptotic Analysis of Population Dynamics
As previously mentioned, we assume that the values of the time intervals {Δ −1 ,  ∈ N}, which the predator spends on search for prey, have an exponential distribution with a parameter  −1 ; that is, E{Δ −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 {  (),   (),  ≥ 0} given by the infinitesimal operator (12).At first it is necessary to check whether such a limit operator  0 exists, that, for any continuously differentiable function V(, ) in some vicinity  0 of each point from the positive quadrant { ≥ 0,  ≥ 0}, the following equality is true: In our case, this limit exists and it is equal to Operator ( 14) corresponds to a deterministic dynamical system of differential equations coinciding with the system of (1).Consequently [22], it can be asserted that for any positive number  there exist such positive numbers  0 and  that for all  <  0 the following inequality exists: sup 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 (+) < (−), 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 {  (0),   (0)} for the system of ( 4)-( 5)-( 6) are chosen sufficiently close to the coordinates {, } of the aforementioned cycle, then the line {  (),   (),  ≤  −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 {  (),   ()} from the mean values {(), ()} for each { ∈ (0,  −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: By definition (17) a two-dimensional random process {  (),   ()} on each interval  ∈ ( −1 ,   ),  ∈ N, is given by the equations but at the moments of contacts   ,  ∈ N, between a predator and a prey it is given by the condition of a jump Substituting in formulas ( 18)-( 19) the expression {  (),   ()} with the expression we can obtain the following system of equations: The process {(), (),   (),   ()} depends only on its state at the time moment ; that is, by definition (17) the two-dimensional random process {  (),   ()} is filtration {F  ,  ≥ 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 {(), (),   (),   ()}.A weak infinitesimal operator L() [18] of this process can be found by calculating the limit expression for an arbitrary finite sufficiently smooth function (, , , ).Considering the possibility of contact between populations in a sufficiently short time interval (0, ), the following expressions can be used to calculate the limit ( 22): Besides, Substituting the asymptotic by  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 {(), (),   (),   ()}: 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 (, , , ) and further use the asymptotic equalities: Using these expansions, formula (25) can be rewritten in the form convenient for passage to the limit: It follows that the limit operator L = lim →0 L() corresponding to the Markov process {(), (), (), (),  ≥ 0}, given by a system of ordinary differential equations (15), is a stochastic linear inhomogeneous stochastic differential equation and an ordinary homogeneous differential equation where and () 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 {(), (),  ≥ 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 {  (),   ()} of the impulse-differential equations system (4)-( 5)-( 6) with accuracy of the order of smallness greater than √ can be approximated by random processes {() + √(), () + √()}.In Figure 3 a sample trajectory of the process {() + √(), () + √(),  ∈ [0, 20]} is simulated in MATLAB environment with the initial conditions (0) = 20, (0) = 3, (0) = 0, (0) = 0, var{(0)} = 0, var{(0)} = 0, covar{(0), (0)} = 0, and  = 0.01 on the phase plane {, }.
Comparison between Figures 1 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 let us rewrite the system of ( 28)-(29) in a vector-matrix form: 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 where () is the matrix solution of the ordinary differential equation with respect to the initial value (0) = ( 1 0 0 1 ).By the definition of a stochastic integral [18], formula (34) defines a Gaussian process with zero mathematical expectation and a covariance matrix The with the initial condition (0) = 0. Figures 4 and 5 show the time evolution of the variances  11 () = E{ 2 ()} and  22 () = E{ 2 ()}.

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 twodimensional Gaussian Markov process.This process has an unboundedly increasing variance, leading to possible significant deviations even on relatively short time intervals.