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.
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]: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 , 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 . 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 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].
2. The Model
Let us propose a model where the contacts between predator and prey occur very often at random time moments , which are the moments of discontinuity of the trajectories of a stationary piecewise constant Poisson process [18] given on a certain probability space with an exponentially distributed length of the constancy intervals , where , and is a small positive parameter. Let us denote as the minimal filtration [19] to which the process is adapted. Let us also assume that at time moments the random variables have a continuous distribution and, therefore, without any loss of generality, we may consider that this distribution is uniform on the interval . It is known [18] that the probability distributions of a homogeneous Markov process are uniquely determined by its weak infinitesimal operator given byIn our case the Markov process is given by an infinitesimal operatorwhere 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 , and densities of prey and predator populations at this moment are , respectively. Suppose that the size of the prey’s population changes in accordance with the logistic equationAnd the size of the predator’s population changes by the rulewhere , , , and . At the time moment the collision between predator and prey individuals occurs, and predator’s gain is given by the expression Here and further . Up to the time moment , the dynamics of the size of the prey and predator populations are also given by (5)-(6) with initial conditions and . At the moment the predator finds the prey again and its gain is given by expression (6) where the argument is replaced by . Then up to the moment 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 equationswhere is a Poisson random measure with the parameter .
Figure 1 shows the trajectory of the solution of the system of the impulse-differential equations (4)-(5)-(6) with the initial conditions and and the parameter values , , , , , , and . Figure 2 shows the trajectory of the solutions of the system of (1) with initial conditions and and with the same the parameter values , , , and .
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 is filtration adapted, and for any and we have the identity:i.e., has the Markov property [18]. In addition, the dynamic characteristics (4)-(5)-(6) of the process are homogeneous in time. Consequently, the dynamical system defined by the Poisson process and (4)-(5)-(6) defines a homogeneous two-dimensional Markov process with the infinitesimal operator [18] given byIn a sufficiently small time interval the stationary Poisson process can perform jumps from point to point , and both random variables are independent and equally uniformly distributed on the segment . Therefore, when the following asymptotic equalities hold:where is an infinitely small quantity with order higher when . Consequently, if is a finite sufficiently smooth function, thensubstituting these asymptotic equations into formula (9), one can obtain an expression for a weak infinitesimal operator for a fixed :
3. Asymptotic Analysis of Population Dynamics
As previously mentioned, we assume that the values of the time intervals , which the predator spends on search for prey, have an exponential distribution with a parameter ; that is, . 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 given by the infinitesimal operator (12). At first it is necessary to check whether such a limit operator exists, that, for any continuously differentiable function in some vicinity of each point from the positive quadrant , the following equality is true:In our case, this limit exists and it is equal toOperator (14) corresponds to a deterministic dynamical system of differential equationscoinciding with the system of (1). Consequently [22], it can be asserted that for any positive number there exist such positive numbers and that for all the following inequality exists: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 for the system of (4)-(5)-(6) are chosen sufficiently close to the coordinates of the aforementioned cycle, then the line , 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 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 , is given by the equationsbut at the moments of contacts , between a predator and a prey it is given by the condition of a jumpSubstituting in formulas (18)-(19) the expression with the expressionwe 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 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 [18] of this process can be found by calculating the limit expressionfor an arbitrary finite sufficiently smooth function . Considering the possibility of contact between populations in a sufficiently short time interval , 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 corresponding to the Markov process , given by a system of ordinary differential equations (15), is a stochastic linear inhomogeneous stochastic differential equationand an ordinary homogeneous differential equationwhereand 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 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 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 is simulated in MATLAB environment with the initial conditions , , , , , , , and 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 . This means that at each time moment the centered second moments of the process can be approximated using the second moments of the process in the form of asymptotic equalities:where , and . For further analysis, using the notationslet 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 integralwhere is the matrix solution of the ordinary differential equationwith respect to the initial value . By the definition of a stochastic integral [18], formula (34) defines a Gaussian process with zero mathematical expectation and a covariance matrixThe covariance matrix that corresponds to the integral (36) as a function from parameter can be represented in the form of integralor in the form of a matrix solution of a system of ordinary differential equationswith the initial condition . Figures 4 and 5 show the time evolution of the variances and .
Figure 6 shows the time evolution of the covariance .
Figures 4–6 show that the variances of the normalized deviations increase very rapidly and even at it cannot be recommended to use the most popular applied statistics “3 rule”; that is, the inequality because the quantities either or can become negative, which contradicts the definition of random processes and . For clarity, Figure 7 shows the process trajectory on the time interval on the phase plane .
Here the variable at some time point already exceeds the value , although the value (the abscissa of the limit cycle in Figure 2) does not exceed 24.
Unbounded increase of variances is a consequence of a resonance in the equation (38). Let be a period of the limit cycle . Then and are -periodic matrix functions and their elements may be decomposed into a Fourier series by and , . Therefore, the integrand in (37) contains a constant component and . At a finite time interval that depends on we can use the diffusion approximation for estimation of finite dimensional distributions of the process .
Even if a deterministic limit cycle closely approaches one of the axes (Figure 8) we can use diffusion approximation (Figure 9) for prognoses of the population growth up to time .
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.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.