Stochastic Analysis of Predator–Prey Models under Combined Gaussian and Poisson White Noise via Stochastic Averaging Method

In the present paper, the statistical responses of two-special prey–predator type ecosystem models excited by combined Gaussian and Poisson white noise are investigated by generalizing the stochastic averaging method. First, we unify the deterministic models for the two cases where preys are abundant and the predator population is large, respectively. Then, under some natural assumptions of small perturbations and system parameters, the stochastic models are introduced. The stochastic averaging method is generalized to compute the statistical responses described by stationary probability density functions (PDFs) and moments for population densities in the ecosystems using a perturbation technique. Based on these statistical responses, the effects of ecosystem parameters and the noise parameters on the stationary PDFs and moments are discussed. Additionally, we also calculate the Gaussian approximate solution to illustrate the effectiveness of the perturbation results. The results show that the larger the mean arrival rate, the smaller the difference between the perturbation solution and Gaussian approximation solution. In addition, direct Monte Carlo simulation is performed to validate the above results.


Introduction
The Lotka-Volterra (LV) model [1] is a classic model for describing the interaction between the prey and predator in ecosystems with predator-prey relations [2,3]. Since the introduction of this model, several improved LV models have been investigated with consideration of different factors. Among those studies, the ones considering predator saturation and predator competition terms, respectively, have been introduced to characterize two special cases of the evolution of the ecosystems. The model with the predator saturation term describes a case where the preys are abundant, while the one with the predator competition term considers another case in which the predator population is very large compared with the prey population. The dynamics of these two models have been discussed in [1,4].
Moreover, the evolution of real-world biological systems usually suffers from unavoidable random perturbations in the natural environment [5][6][7][8][9]. Many investigations on the dynamics of predator-prey models with predator saturation and predator competition terms excited by stochastic perturbations have been reported [8][9][10][11][12][13][14]. In these studies, the stochastic perturbation is introduced as continuous stochastic excitations to disturb the birth rate of the prey and the death rate of the predator.
Nevertheless, in the environmental fluctuations, there are some drastic changes, including floods, earthquakes, tornados, or forest fires, which indicates that the evolutions of the ecosystem can be affected by unavoidable, sparse, discrete random jumps. Consequently, it is not appropriate to characterize the environmental effects by using only continuous random excitations. The combinations of continuous random processes and random jumps are viewed as more appropriate model for this case [15][16][17][18]. The ecosystem models excited by stochastic excitations with random jumps have been attracting increasing interests, and many reliable results have been obtained. In these studies, the Lévy process or jump-diffusion processes are commonly used in mathematical models to describe environmental perturbations with random jumps [19][20][21][22]. Most of them focus on the mathematical properties of the ecosystems, including the uniqueness, stability of the solutions, etc. [23][24][25][26][27]. However, for the prediction of the PDFs of the species populations, only the ecosystems enforced by stochastic processes with jumps have been investigated [28,29]. Therefore, more work needs to be done on deriving the stochastic responses of the ecosystems under stochastic excitations with random jumps.
The stochastic averaging method is an effective tool for investigating nonlinear systems excited by stochastic forcings. It was first applied to nonlinear systems under Gaussian white noise excitation [30][31][32], and then generalized to nonlinear systems with other types of stochastic excitations [33,34]. It has been successfully employed to study the stationary PDFs [6][7][8][9]11,12,28,35] and the optimal control [36] of species populations in ecosystems with small self-competition and small stochastic excitation. However, the solution of the ecosystem dynamics under continuous and random jump excitation by the stochastic average method has not yet been reported.
Inspired by this, here, we study predator-prey models with predator saturation and predator competition terms excited by combined Gaussian and Poisson noises. The governing equation for the solution of the stochastic model is approximately derived by applying the stochastic averaging method. Then, the statistical responses including the approximate PDFs and moments of the population densities are derived by solving the governing equation with a perturbation technique. Further, the effects of the system parameters as well as the stochastic excitations on the stationary responses are presented. The rest of the paper is organized as follows. Section 2 introduces two ecosystem models excited by combined Gaussian white noise and Poisson white noise and the averaged generalized Fokker-Planck-Kolmogorov (GFPK) equation. In Section 3, the effects of the parameters on the stochastic dynamics of the ecosystem are discussed. Finally, the conclusions are summarized in Section 4.

The Deterministic Models
In the present paper, we start from two deterministic prey-predator type ecosystems, one with an abundant prey population and the other with a large predator population, which can be viewed as generalizations of the LV model [1]. Although the two models are given based on different assumptions, fortunately we can define some new parameters to obtain a more general model that contains the two models as its special cases [8,37].

Case 1: prey population is abundant
Consider the case where the prey is abundant. The model is given as .
where x 1 and x 2 are the prey and predator population densities, respectively; a and c represent the birth rate of the prey and the death rate of the predators. −sx 2 1 reveals the self-competition between the preys. The terms bx 1 x 2 1+Ax 1 and f x 1 x 2 1+Ax 1 in Equation (1) describe the interaction between the species [37]. These two interaction terms are proportional only to the predator population x 2 , which implies that the consumption of the prey depends on the predator population and not on the prey population.
Model (1) can be reformulated in the following form . with by introducing At this time, Model (2) has two equilibrium points, one unstable point (0, 0) and one asymptotically stable point c/ f , a/b .

Case 2: the predator population is large
When the population predator is large relative to the prey population, the preypredator system can be modeled as .
In this case, the prey consumption depends only on the prey supplement. Compared with model (1) in Case 1, Model (5) has different interaction terms bx 1 x 2 1+Bx 2 , f x 1 x 2 1+Bx 2 , which imply that the prey growth rate and the predator death rate depend only on the prey supply. [8]. The other terms in (5) are the same as those in Equation (1).
By defining the new parameters Model (5) can be expressed in the form of Equation (2) with different nonlinear functions: Some assumptions about the systems in Case 1 and Case 2 are needed for subsequent analysis. First, we consider the situation that the parameter s in the self-competition term for both cases is small. It describes the situation that the prey population density is small, and the impact of the self-competitions term −sx 2 1 is small. Second, we assume the term Ax 1 1 for the Case 1. This assumption is consistent with the first one. Third, for Case 2, Bx 2 1 is also necessarily needed to meet the requirement of the first assumption. For the convenience of following analysis, a small and positive parameter ε is introduced and the parameters s, A and B in Equation (2) are replaced by ε 2 s, ε 2 A, and ε 2 B to represent small values of these parameters. Thus, the mathematical model that we will investigate becomes . with for Case 1 and with for Case 2.

The Stochastic Model
We take the random environmental fluctuations into account in Equation (8). Gaussian white noise and Poisson white noise are adopted to model the random environmental fluctuations, which cause random changes in the growth rate of the prey and death rate of the predators. This means that the growth rate of the prey and the death rate of the predators in Equation (8) are changed as: Here, ε is the small parameter from Equation (8), which indicates that the random environmental influences are small. ζ i (t)(i = 1, 2) are two independent Gaussian white noise excitations that are used to characterize the continuous environmental fluctuations. They have the following statistical characteristics: where E[•] represents the mathematical expectation and 2D i are the noise intensities. Moreover, two independent Poisson white noises ξ i (t)(i = 1, 2) are used to model the jump effects in the environment, which can be viewed as the formal derivatives of compound Poisson processes C i (t) In Equation (14), N i (t) denote Poisson counting processes with mean arrival rates λ i > 0. Y ik represent the random magnitudes of the impulses that obey Gaussian distributions with variances E[Y 2 ik ] and zero mean in the present paper. U(•) is the unit step function.
ik ] are the intensities of the Poisson white noises. In addition, ζ i (t) are also assumed to be independent with ξ i (t).
Consequently, the system that we will investigate can be described by the following Stratonovich stochastic differential equation (SDE) based on Equations (8) and (11), where for Case 1 and for Case 2. Mathematically, the symbols X 1 (t) and X 2 (t) are used to represent the stochastic processes corresponding to the deterministic functions x 1 (t) and x 2 (t) in Equation (8).
Equation (14) can be transferred to the following Itô SDE by adding some correction terms as follows [38] in which X 1 i! (dC 2 (t)) i are two correction terms. We introduce the integral form of the compound Poisson process by Poisson random measures for further analysis [39] where P k (dt, dY k ) are Poisson random measures and Q k are the Poisson mark spaces. Based on the Poisson random measure, the differentiation of compound Poisson processes can be written as Then, the equivalent stochastic integro-differential equation (SIDE) of Equation (17) can be derived as [40] where In the remainder of this paper, the stochastic dynamics of the stochastic ecosystems will be analyzed based on Equation (20). Different from the deterministic system, the probability measures for the prey and predator population are needed to describe the random behavior of the system with random excitations. For the study of the stochastic problem, probability plays an important role. Thus, an important task here is to obtain the PDFs for both prey and predator population densities. In this paper, we focus on the longterm behavior of this stochastic model, namely, the stationary responses of Equation (14). To derive the stationary PDFs, we use a stochastic averaging method as described below.

The Stochastic Averaging
Without considering the random excitations and the nonlinear function g i (x 1 , x 2 ) (I = 1, 2) and letting the self-competition term be 0, Equation (20) reduces to the famous LV model .
It is known [1] that the system (22) has periodic trajectories that can be defined by where Q is a constant and r(x 1 , x 2 ) = 0 at both the unstable equilibrium point (0, 0) and the stable non-asymptotic equilibrium point c/ f , a/b of Equation (22). This means that Equation (23) is the first integral of the system (22) for any positive value Q [1]. The period of the periodic trajectories with constant Q is determined by Figure 1 demonstrates the trajectories of the system for three cases. Figure 1a shows the periodic trajectories for system (22) with different initial values. For different initial values, different closed orbits are obtained, namely, periodic solutions. Figure 1b,c show the random trajectories for system (22) under Gaussian white noise and Poisson white noise, respectively. It is known that stochastic excitation has a significant influence on the dynamics of the ecosystems. The trajectories of the system under stochastic excitations are random. Moreover, the trajectory changes with continuous slight jumps for the system with Gaussian white noises, while for the system excited by Poisson white noises, the trajectory includes some large jumps, which demonstrates the influences of different stochastic excitations on the trajectories. It can also be seen in these figures that large population densities of species may evolve into smaller population densities with the evolution of time. For an ecosystem, a smaller population density means that the population is more likely to die out, which means that the system is relatively fragile or unstable. Replacing the variables x 1 and x 2 by stochastic processes X 1 and X 2 in Equation (23), the random counterpart of the first integral (23) can be obtained as Based on the stochastic chain rule [39], the differential form of the first integral is given as Taking the Taylor expansions of ln(1 + γ 11 /X 1 ) and ln(1 + γ 22 /X 2 ), and substituting them in Equation (26), one arrives at in which 24 . Replacing X 1 in Equation (20) by R(t) yields a new system governed by Equation (27) and the second equation of Equation (20). R(t) is a slowly varying variable since ε is a small parameter, while X i are rapidly varying variables. Based on the stochastic averaging method proposed by [41,42], one can derive the averaged GFPK equation for R(t) as where p(r, t) is the PDF of R(t) at time t. O(ε 5 ) represents the terms of ε 5 and higher, which contains an infinite number of terms of GFPK equation with small magnitude due to the increasing powers of ε. The other parameters of Equation (28) are given as where the symbol [ ] t means the time average in one quasi-period, which is defined as where T is the period given in Equation (24). It can be found from Equations (29)-(37) that the coefficients consist of terms of order ε 2 and ε 4 . In the absence of terms of order ε 4 and higher in Equation (28) the averaged GFPK Equation (28) will reduce to the averaged FPK equation for the system under only Gaussian white noises ζ i (t) with intensities (2D i + λ i E[Y 2 i ])(i = 1, 2). By solving this reduced FPK equation, one can obtain the Gaussian approximation solution for the averaged GFPK equations. However, the Gaussian approximation solution is not a proper approximation for small values of λ for Poisson white noise with the same noise intensity since the influence of term of order ε 4 in Equation (28) cannot be ignored. Therefore, more terms of the GFPK equations should be considered to obtain a more accurate solution when dealing with the system with Poisson white noise.

The Approximate Stationary Responses
The PDFs of the population densities are obtained by solving the averaged GFPK Equation (28) with certain boundary and initial conditions. In the present paper, only the long-term behaviors of the ecosystem are studied. A perturbation technique [43] is applied to derive the stationary PDFs p R (r) for R(t) by solving the following reduced averaged GFPK equation A second-order perturbation solution p R (r) = p 0 (r) + εp 1 (r) + ε 2 p 2 (r) (40) is adopted here to derive the approximate solution of Equation (39). By substituting Equation (40) into Equation (39) and putting terms of the same order of ε together, the equations for p 0 (r), p 1 (r), p 2 (r) are given as It can be seen from Equation (41) that p 0 (r) is the solution for the system excited by Gaussian white noise with the same noise intensity, namely, approximate Gaussian solutions for the system excited by Poisson white noise. By solving Equations (41)-(43) step-by-step, the second-order perturbation solution (40) can be obtained. Then, the approximate stationary joint PDFs p X 1 X 2 (x 1 , x 2 ) take the form The marginal PDFs and moments of X 1 and X 2 can be calculated from Equation (44), such as Furthermore, when the Poisson white noise intensity is 0, the method developed in this paper reduces to the case excited by Gaussian white noise. The effects of the Gaussian noise on the dynamics of the ecosystems have been studied by Cai [7]. Therefore, in the following subsections, we focus on the influence of the system parameters and Poisson white noise parameters on the stationary statistics.

The Influence of System Parameters
In this section, the influences of system parameters on the model (14) for Case 1 and Case 2 are shown in Figures 2-9. The stochastic properties of population densities of the prey and predator can be calculated from Equations (44)-(47).        Case 1: prey is abundant compared with the predator The effects of the parameter ε 2 A in the predator saturation term of Equation (15) on the stationary PDFs of the species are shown in Figure 2. The PDF p X 1 (x 1 ) of the population density of the prey X 1 (t) and the PDF of the population density of the predators X 2 are calculated with the following parameters a = 1.
i ] = 0.02, i = 1, 2, and two different values of ε 2 A. One can see from Figure 2 that the occurrence of the predator saturation influences the system behavior significantly. Compared with the case of ε 2 A = 0.0, the maxima of the PDFs of the species population densities become lower, and the probabilities for both small and large populations increase for ε 2 A = 0.05. This means that the ecosystem becomes more unstable for ε 2 A = 0.05. This is reasonable in the real world. Since the parameter ε 2 A represents the case of large prey population, it is obvious that the ecosystem is more unstable when this extreme case occurs. The solid lines in Figure 2 show the second-order perturbation solutions, while the discrete points denote the results obtained from Monte Carlo simulation. The good agreement between these two results shows the effectiveness of the second-order perturbation solution.
The effect of ε 2 A on the moments of the population densities of prey and predators which can be calculated from Equations (46) and (47) are depicted in Figure 3. It can be seen that the curves for the means and relative fluctuations Var(X i )/E[X i ] of the population densities of both species increase monotonously. This means that the ecosystem fluctuates more strongly for larger value of ε 2 A, which implies a more unstable system. Case 2: predator population is large Figure 4 depicts the PDFs of X 1 and X 2 for the model in Case 2 for different ε 2 B values. It is found that the maxima of the PDFs for the case of ε 2 B = 0.05 are larger and the probability in this case is more concentrated. Figure 5 shows the mean and relative fluctuation of the species population of the ecosystem. It can be seen that the mean value of the prey and predator increase with increase of the ε 2 B value, while the relative fluctuation curves decrease with increasing the value of ε 2 B, implying a more stable ecosystem. The results of the Monte Carlo simulations are also given in Figures 4 and 5 to validate the results of the proposed method.

The Influence of Poisson White Noise
In this section, we focus on the influences of the parameters of the Poisson white noises, namely mean arrival rate λ and the variance of amplitude E[Y 2 ].
In Figures 6 and 7, some numerical calculations have been performed to obtain the PDFs of predator and prey population densities in the model (14) for Case 1 and Case 2, respectively. In these figures, the second-order perturbation results are represented by the solid lines. The dashed lines are the approximate Gaussian solutions. The Monte Carlo results are plotted as dotted lines. It is found that the PDFs obtained by the present method are closer to the Monte Carlo simulation than the approximate Gaussian solutions. The maxima of the stationary PDFs for the system with both Gaussian and Poisson white noise are higher than for the approximate Gaussian solutions. This means that the second-order perturbation solutions have higher accuracy than the Gaussian approximation solution. Figures 8 and 9 show the influences of λ of the Poisson white noises on the stationary PDFs when the Poisson white noise intensity λE[Y 2 ] is kept constant. It is observed from these figures that with the increase of λ, the maxima of the PDFs decrease and the PDFs for population densities of prey and predator approach the ones for the system with the same intensity of Gaussian white noise excitation. To see this clearly, we define the L 2 error norm between the PDFs for the system with Poisson white noise and the Gaussian approximation solution as where p X i (x i ) are the second-order perturbation PDFs for X i and p G x i (x i ) are the Gaussian approximation solutions. We calculate the errors for different λ for Case 1 ( Table 1) and Case 2 ( Table 2). It can also be seen from these tables that with the increase of λ, the error decreases significantly, which also verifies our above conclusion.

Conclusions
In the present paper, we have investigated the statistical responses of a stochastic prey-predator type model for the possible situations of sufficient prey supply and a large predator population. The statistical responses of the species population, including the approximate stationarity PDFs and moments, have been obtained by a stochastic averaging and a perturbation technique. It can be found that, for the system with abundant prey, the increase of the parameter ε 2 A describing the saturation of the predator population from 0 to 0.5 leads to the increase of population fluctuations, which implies that the system becomes unstable. However, for the system where the predator population is large, the increase of the parameter ε 2 B describing the competition of predators for prey has the opposite effect. In addition, we have paid special attention to the effect of the mean arrival rate λ of the Poisson white noise. The influence of Poisson white noise on the system tends to approach the influence of Gaussian white noise with the same noise intensity when the λ increases from 0.05 to 2.0.
Although only two types of ecosystems with combined Gaussian and Poisson white noise have been investigated in this paper, our approaches can be applied to other ecosystem models. Except for the stochastic response of the ecosystem being investigated, the stochastic optimal control or extinction problem are also worth studying.