Stochastic Dynamics of a Time-Delayed Ecosystem Driven by Poisson White Noise Excitation

We investigate the stochastic dynamics of a prey-predator type ecosystem with time delay and the discrete random environmental fluctuations. In this model, the delay effect is represented by a time delay parameter and the effect of the environmental randomness is modeled as Poisson white noise. The stochastic averaging method and the perturbation method are applied to calculate the approximate stationary probability density functions for both predator and prey populations. The influences of system parameters and the Poisson white noises are investigated in detail based on the approximate stationary probability density functions. It is found that, increasing time delay parameter as well as the mean arrival rate and the variance of the amplitude of the Poisson white noise will enhance the fluctuations of the prey and predator population. While the larger value of self-competition parameter will reduce the fluctuation of the system. Furthermore, the results from Monte Carlo simulation are also obtained to show the effectiveness of the results from averaging method.


Introduction
In the real-world ecosystem, a time delay effect in the ecosystems is inevitable. Such as, in the predator-prey type ecosystem, it takes time for the predator population to adjust any change of the prey population. Thus, the ecosystems with time delay effects have received more and more attentions in past decades. Several types of models considering the time delay effects are introduced. The dynamical properties of the ecosystem with time delay, such as stability and bifurcations, have been studied [1][2][3][4][5][6][7]. It indeed can be found that the effect of the time delay can change the behavior of the ecosystem substantially. However, the models in all these studies are deterministic. They are too idealistic since the changes in the environment are not taken into consideration. As a matter of fact, random perturbation is universal in the real world. Thus, the stochastic excitations describing the fluctuations in the environment are added to the deterministic models. The dynamics of the ecosystem under stochastic excitations have been widely investigated [8][9][10][11][12][13][14][15][16][17][18][19][20]. Especially, the stochastic dynamics of the noise-induced predator-prey type ecosystems with time delay, including stochastic responses [9,19,21,22], stochastic stability [23], stochastic resonance [16,24] and noise-delayed extinction [16], have been investigated by some numerical and analytical methods.
It can be noted that the random perturbations in the environment in the previously mentioned works are usually modelled as continuous stochastic processes, such as the Gaussian white noise or Gaussian colored noise. However, the perturbations in the real-world environment are complicated. There are some unavoidable sparse, drastic changes in the environment, for instance sudden natural disasters, earthquakes, forest fires, floods [25][26][27]. These changes are pulse-type perturbations, which cannot be characterized by the continuous stochastic processes properly. Therefore, the stochastic jump The second equation of Equation (1) can also be written as following equivalent form [3] .
where x 1 and x 2 are the population densities of prey and predator, respectively. The parameters a 1 , s, b, c and f are positive constants. a 1 is the birth rate of the preys. c is the death rate of the predators. The term −sx 2 1 denotes the prey self-competition. The terms −bx 1 x 2 and f x 2 t −∞ F(t − τ)x 1 (τ)dτ provide a balance between the prey and predator populations. The integral term represents the time delay effect of the prey population on the predator population. It means that the change of the prey population affects the predator population after a time lag. And this time delay depends on an average over past populations, not only on the population at some particular instant in the past. F(·) in the integrand is called the delay function and normalized as follows:  Macdonald has presented two ways to reduce the system (1) in [3]. The first one is exact, but, it depends on the choice of a simple form of the time delay function F(·). The other one does not need to specify the time delay function, but needs its moments, i.e.: Equation (3) can be viewed as a measure of the average delay time. Two reasonable choices for F(t) provieded by Cai and Lin [9] are: γ e −t/γ and F(t) = 4 To avoid specifying time delay function, the second way is used to reduce the system (1) in the present paper. As pointed out in [3], this method is valid for small γ. Substituting the first-order approximation of the Taylor expansion of x 1 (t − τ) about τ = 0 into Equation (2) and applying Equation (3), the first-order approximation of Equation (2) can be derived as: Together with the first equation of Equation (1), the original Equation (1) can be rewritten as: .
Let a = a 1 − sc/ f . Then the system can be converted to: .
The system (6) can be used to investigate the effect of the time delay on the populations of the prey and predator. It is easy to check that system (6) has the same equilibrium point x 10 = c/ f , x 20 = a/b as that of the Lotka-Volterra model without time-delay and stochastic excitations [8,9]. In addition, it is pointed in Cai and Lin's work [9] that, for a physically meaningful and stable ecological system, the choice of the parameters must make the ecosystem stable. All the parameters in the present paper meet this requirement.

Stochastic Model
In the present paper, we suppose that there are some unavoidable sparse, drastic changes in the environment. These changes are modeled as the Poisson white noises. They will cause random jump in the prey growth rate and the predator death rate. Thus, the stochastic model has following form: where W 1 (t) and W 2 (t) are two independent Poisson white noises, which can be defined as [35,36] where the δ(•) is the Dirac delta function; N i (t) is a Poisson counting process with mean arrival rate λ i > 0 and gives the number of the pulses arriving in the time interval [0, t]; Y ik represents the amplitude of the Poisson white noises, which is independent of the pulse arrival time t ik . In the present paper, Y ik is Gaussian distribution with mean value 0 and variance E[Y 2 i ] [32,33]. It is hard to obtain the exact solution of Equation (7) analytically, even for the Gaussian white noise case. Therefore, some methods have been developed to solve this problem. Among which, the stochastic averaging method is an analytical efficient method. It has been successfully applied to analyze the stochastic behaviors of the ecosystem under stochastic continuous excitations [18,19,22,37]. In recent years, the stochastic averaging method has been generalized to investigate the nonlinear system with random jump excitations [38][39][40][41]. In the following section, we will employ the stochastic averaging method to study stochastic dynamics of the ecosystem (7).
To proceed further analysis, the following assumptions about the system are in order. First, the coefficient s in self-competition term is generally small. It means that the self-competition term has small influence when the prey population density is small. Second, the parameter γ is small, corresponding to shorter time-delay effect. Finally, the noise intensities of W 1 (t) and W 2 (t) are small, indicating small random fluctuations in the system dynamics. These assumptions are generally valid for the real-world ecosystem [8,9]. For the convenience of the further theoretical analysis, in the following part, the parameter s, γ and λ i E[Y 2 i ] in Equation (7) are replaced by: ε 2 s, ε 2 γ and λ i ε 2 E[Y 2 i ], respectively, where ε is a small parameter. Thus, Equation (7) can be rewritten as: Equation (9) is usually modeled as the Stratonovich stochastic differential equation (SDE) [9,13,18]. By adding some correction terms [13,35,36], it can be converted to following equivalent Itô SDE: The terms X 1 (10) are the correction terms from Stratonovich SDE to Itô SDE. Then, Equation (10) can be converted to following equivalent stochastic integro-differential equation (SIDE) [42]: in which: and P 1 (dt, dY 1 ) (i = 1, 2) are the Poisson random measures; Y i (i = 1, 2) are Poisson mark spaces [42].

Stochastic Averaging Approach
In this subsection, a stochastic averaging procedure is derived to simplify the original system. For future use, we choose the following deterministic conservative system for the stochastic averaging [9,18]: Entropy 2018, 20, 143 5 of 15 The system (13) has the same equilibrium point as that of system (9) without the Poisson white noises. System (13) possesses a first integral [9] r(x 1 , This first integral can be view as an extension of first integral for the standard LV system [8,43]. It can be seen that r(x 1 , x 2 ) = 0 at the equilibrium point (c/ f , a/b). And for each positive constant K, r(x 1 , x 2 ) = K represents a periodic trajectory surrounding the equilibrium point. For illustrative purpose, Figure 1 shows the equilibrium point O and three different periodic trajectories for system (13) with ε = 0.1, a = 0.9, b = 1.0, c = 0.5, f = 0.5, ε 2 s = 0.2, ε 2 γ = 0.2. In this figure, the equilibrium point corresponds to K = 0, and the trajectories correspond to K = 0.892, 1.146, 1.565. And the period of the trajectory can be determined from in which the integration symbol means integration along the trajectory r(x 1 , The system (13) has the same equilibrium point as that of system (9) without the Poisson white noises. System (13) possesses a first integral [9]     (14) This first integral can be view as an extension of first integral for the standard LV system [8,43].
It can be seen that   , r x x K represents a periodic trajectory surrounding the equilibrium point. For illustrative purpose, Figure 1 shows the equilibrium point O and three different periodic trajectories for system . In this figure, the equilibrium point corresponds to  0 K , and the trajectories correspond to  ， 0.892 K 1.146 , 1.565 . And the period of the trajectory can be determined from x a bx x c fx ε sγx (15) in which the integration symbol means integration along the trajectory    1 2 , r x x K. Replace 1 x and 2 x by the stochastic processes   1 X t and   2 X t , the stochastic counterpart of can be given as: Based on the stochastic jump-diffusion chain rule [42] and Equation (11), the SIDE for   R t is: Replace x 1 and x 2 by the stochastic processes X 1 (t) and X 2 (t), the stochastic counterpart of Equation (14) R(t) = R(X 1 (t), X 2 (t)) can be given as: Based on the stochastic jump-diffusion chain rule [42] and Equation (11), the SIDE for R(t) is: Entropy 2018, 20, 143 6 of 15 In order to possess the stochastic averaging method, consider the Taylor expansion of ln(1 + h 11 /X 1 ) and ln(1 + h 22 /X 2 ) in Equation (17), and then substitute Equation (12) to Equation (17). After collecting the terms of same order of ε, Equation (17) can be written as: where: It can be seen from Equation (18) that R(t) is a slowly varying stochastic process since ε is a small parameter. According to the stochastic averaging method for the nonlinear system under the Poisson white noise [33], the averaged GFPK equation can be derived as: where O ε 6 contains the term of the order of ε 6 and higher, and the coefficients for the GFPK equation are given in Appendix A. Note that the SIDE (18) contains infinite terms. To get the closed form of averaged equations, truncation is needed during the averaging process. In the averaged GFPK Equation (25), only the terms up to the fourth order of ε are considered for simplicity.

Stationary Probability Density Functions
The averaged GFPK equation can be solved with certain boundary and initial conditions. The stationary PDF of R(t), denoted by p(r), can be obtained by solving the reduced GFPK equation: with the following conditions: The perturbation method [33,44] is used to solve the reduced GFPK equation. The second order perturbation solution p(r) = p 0 (r) + εp 1 (r) + ε 2 p 2 (r) is adopted in our calculation. On substituting Entropy 2018, 20, 143 7 of 15 the second order perturbation solution into Equation (27) and grouping terms of the same order of ε, the following set of differential equations is obtained: In In ε 4 , Equation (28) is a Fokker-Planck-Kolmogorov equation for Gaussian white noise excitation. The exact solution can be obtained if the FPK Equation (28) belongs to the class of generalized stationary potential [39]. The perturbation solutions p 0 (r), p 1 (r) and p 2 (r) can be obtained by solving Equations (28)-(30) successively.
After getting p(r), the stationary PDFs of X 1 and X 2 can be calculated as follows [9,13]: where r is the function of x 1 and x 2 given in Equation (14) and T(r) is the quasi-period with the form in Equation (15). The other statics can be easily derived from Equation (31).

Results
Some numerical results are presented in Figures 2 and 3 for ecosystem (7) with different noise intensities. The parameters are shown in the captions of these figures. In these figures, the solid lines are the results of ecosystem under Poisson white noises obtained by the proposed method. The dotted lines are the results from Monte Carlo simulation. A good agreement between the theoretical and simulation results in both figures shows the accuracy of the proposed method. For the purpose of comparison, the stationary PDFs for the time-delayed ecosystem under Gaussian white noises with the same intensity are also depicted in these figures. They are denoted by the dashed lines. One can see that the approximate stationary PDFs for time-delayed ecosystem under Poisson white noise are higher than those for system under Gaussian white noise. It means that, for the same noise intensity, the influence of the Gaussian white noises on the ecosystem is stronger than that of the Poisson white noises.
In the Monte Carlo simulation for ecosystem (7), the Runge-Kutta method for the Poisson white noise excitation proposed by Di Paola and the triangular pulse model of the Poisson white noise excitation are adopted [35]. It is important to note that the numerical simulation time should be long enough to get the stationary responses. In present paper, for each set of system parameters, eight samples are simulated. The time of Monte Carlo simulation is 400,000,000 units for each sample, and the time step length is 0.001.
Cai and Lin [8] show that the stability of system (6) depends on the prey self-competition parameter and the time-delay parameter. Thus, in the following sections, we will investigate the influences of the parameters ε 2 s and ε 2 γ on the ecosystem (7). The stationary PDFs and the statics of prey and predator populations for different values of these two parameters are given and discussed in following sections. In the Monte Carlo simulation for ecosystem (7), the Runge-Kutta method for the Poisson white noise excitation proposed by Di Paola and the triangular pulse model of the Poisson white noise excitation are adopted [35]. It is important to note that the numerical simulation time should be long enough to get the stationary responses. In present paper, for each set of system parameters, eight samples are simulated. The time of Monte Carlo simulation is 400,000,000 units for each sample, and the time step length is 0.001.
Cai and Lin [8] show that the stability of system (6) depends on the prey self-competition parameter and the time-delay parameter. Thus, in the following sections, we will investigate the influences of the parameters 2 ε s and 2 ε γ on the ecosystem (7). The stationary PDFs and the statics excitation are adopted [35]. It is important to note that the numerical simulation time should be long enough to get the stationary responses. In present paper, for each set of system parameters, eight samples are simulated. The time of Monte Carlo simulation is 400,000,000 units for each sample, and the time step length is 0.001.
Cai and Lin [8] show that the stability of system (6) depends on the prey self-competition parameter and the time-delay parameter. Thus, in the following sections, we will investigate the influences of the parameters 2 ε s and 2 ε γ on the ecosystem (7). The stationary PDFs and the statics

The Effects of the Time Delay Parameter ε 2 γ
In this subsection, we investigate the influence of the time delay on the ecosystem (7).

Figures 4 and 5 show the effects of the time delay parameter on the PDFs and relative fluctuations
Var(X i )/E(X i ) of the prey population and predator population. Figure 4a,b depict the stationary PDFs of the prey population p X 1 (x 1 ) and predator population p X 2 (x 2 ) for different values of time delay ε 2 γ, respectively. In these figures, the solid lines are the theoretical results obtained by the proposed method. It can be found that, with increasing the time delay ε 2 γ from 0.05 to 0.12, the peak values of both p X 1 (x 1 ) and p X 2 (x 2 ) become lower, and the probabilities in both lower and higher population become higher. Besides, the results from Monte Carlo simulation are also plotted to show the effectiveness of the proposed method. The other parameters of the system are given as: Figure 5 shows the relative fluctuations of the prey population and predator population changing versus the values of the time delay parameter ε 2 γ for different self-completion parameter ε 2 s. In this figure, the results obtained by the Monte Carlo simulation are represented by point lines, while the results obtained by the proposed method are shown by solid or dashed lines. Each curve of relative fluctuation increases monotonously with the increase of time delay. It means that the fluctuations of the prey population and predator population increase when the time delay is relatively large. Namely, the longer time delay will lead to a wider range distribution of the population, indicating that the ecosystem becomes less stable. of prey and predator populations for different values of these two parameters are given and discussed in following sections.

The Effects of the Time Delay Parameter 2 ε γ
In this subsection, we investigate the influence of the time delay on the ecosystem (7).     of prey and predator populations for different values of these two parameters are given and discussed in following sections.

The Effects of the Time Delay Parameter 2 ε γ
In this subsection, we investigate the influence of the time delay on the ecosystem (7).

The Effects of the Self-Competition Parameter ε 2 s
In this subsection, we will investigate the influence of the self-competition parameter ε 2 s. In Figure 6, the stationary PDFs of prey population and predator population for different values of ε 2 s are given. Figure 7 shows the dependence of the relative fluctuations Var(X i )/E(X i ) on the self-competition parameter ε 2 s for different values of time delay parameter ε 2 γ.
In Figure 6, it can be seen that, for fixed value of other parameters, with the increasing value of ε 2 s, the PDFs become more concentrated, and the peak value of the PDFs become higher. This implies that the system will become more stable when increasing ε 2 s from 0.08 to 0.15. In addition, the results obtained from the Monte Carlo simulation are also calculated to show the validity of the proposed method in the figure. The other system parameters are: ε = 0.1, a = 0.9, b = 1.0, f = 0.5, c = 0.5, Figure 7 shows the relative fluctuations of the prey population and the predator population versus the self-competition parameter ε 2 s for different values of time delay ε 2 γ. Compared with Figure 5, the trends of the relative fluctuations curves are different. For each value of the time delay parameter ε 2 γ, the curve of the relative fluctuations of both prey population and predator population decrease monotonously with the increase of the self-competition parameter ε 2 s. It reveals that the fluctuations of the system will decrease when increase ε 2 s, which indicates that the system becomes more stable for larger ε 2 s. In this subsection, we will investigate the influence of the self-competition parameter 2 ε s . In   In Figure 6, it can be seen that, for fixed value of other parameters, with the increasing value of 2 ε s , the PDFs become more concentrated, and the peak value of the PDFs become higher. This becomes less stable. In this subsection, we will investigate the influence of the self-competition parameter 2 ε s . In   In Figure 6, it can be seen that, for fixed value of other parameters, with the increasing value of 2 ε s , the PDFs become more concentrated, and the peak value of the PDFs become higher. This

The Effects of the Poisson White Noise
It is known that the mean arrival rate λ and the variance of the noise amplitude E[Y 2 ] are two key parameters for the Poisson white noise. In the following section, the influences of λ and E[Y 2 ] are investigated respectively. Figure 8 shows the effects of the mean arrival rate λ 1 = λ 2 = λ on the stationary PDFs. The results are calculated with following parameters: ε = 0.1, a = 0.9, b = 1.0, c = 0.5, f = 0.5, ε 2 s = 0.2, ε 2 γ = 0.1. It can be seen that, when fix the value .01, with increasing the mean arrival rate λ from 0.5 to 2.0, the range of the vibration of the predator population and prey population become larger, and the peak values of the PDFs become lower. And the probabilities in both lower and higher population become higher, indicating a less stable system. This is reasonable. A larger mean arrival rate of the Poisson white noise implies more pulses per unit time. The more the number of pulses is, the more unstable the ecosystem will be. The same conclusion can be obtained from Figure 9, which shows the relative fluctuation of the predator population and prey population versus the mean arrival rate λ. As shown in this figure, the curves of the relative fluctuations for both the predator population and prey population increase monotonously as increasing the mean arrival rate λ. This implies that the system has become unstable. The other parameters are the same as those in Figure 8. Also depicted in Figures 8 and 9 are results obtained from the Monte Carlo simulation.
in both lower and higher population become higher, indicating a less stable system. This is reasonable. A larger mean arrival rate of the Poisson white noise implies more pulses per unit time. The more the number of pulses is, the more unstable the ecosystem will be. The same conclusion can be obtained from Figure 9, which shows the relative fluctuation of the predator population and prey population versus the mean arrival rate λ . As shown in this figure, the curves of the relative fluctuations for both the predator population and prey population increase monotonously as increasing the mean arrival rate λ . This implies that the system has become unstable. The other parameters are the same as those in Figure 8. Also depicted in Figures 8 and 9 are results obtained from the Monte Carlo simulation.  In Figures 10 and 11, the influences of   on the ecosystem are calculated with system parameters: [ ] ε E Y , the fluctuation of the pulses will become larger, which will result in bigger fluctuation of ecosystem. Besides, the results from Monte Carlo simulation represented by dotted lines are also shown in these two figures.
Shown in Figure 12 are the stationary PDFs of the prey and predator populations for different mean arrival rates of the Poisson white noises, corresponding to the same noise intensity  In Figures 10 and 11, the influences of on the ecosystem are calculated with system parameters: ε = 0.1, a = 0.9, b = 1.0, c = 0.5, f = 0.5, ε 2 s = 0.2, ε 2 γ = 0.1, λ 1 = λ 2 = 0.8. Figures 10 and 11 depict the dependence of the stationary PDFs and relative fluctuations of the predator population and prey population. We can arrive at similar conclusions obtained from Figures 8 and 9. As increasing the variance of the impulses ε 2 E[Y 2 ], the ecosystem becomes less stable. This agrees with the intuitive consideration. When increase the value of ε 2 E[Y 2 ], the fluctuation of the pulses will become larger, which will result in bigger fluctuation of ecosystem. Besides, the results from Monte Carlo simulation represented by dotted lines are also shown in these two figures.
Shown in Figure 12 are the stationary PDFs of the prey and predator populations for different mean arrival rates of the Poisson white noises, corresponding to the same noise intensity The results presented in this figure are all obtained by using the proposed method. It is seen that, as increasing the value of the mean arrival rate of the Poisson white noises from 0.3 to 5.0, the approximate stationary PDFs of prey and predator populations of proposed ecosystem under Poisson white noises become closer to those of ecosystem under Gaussian white noises with the same intensity. This implies that the non-Gaussian behavior depends on the mean arrival rate of the Poisson white noise. The other system parameters are given in the caption of this figure.

Conclusions
This paper investigates the stochastic behaviors of a time-delayed predator-prey ecosystem under pulse-type stochastic environment fluctuations. In this model, the stochastic fluctuations are characterized by the Poisson white noises, and the time delay in the interaction between the prey and predator is described approximately by a time delay parameter. The original ecosystem is first modeled as the Stratonovich SDE and then transferred to Itô SDE. Under some reasonable hypothesises, the stochastic averaging method is applied to simplify the original system. Finally, the approximate stationary

Conclusions
This paper investigates the stochastic behaviors of a time-delayed predator-prey ecosystem under pulse-type stochastic environment fluctuations. In this model, the stochastic fluctuations are characterized by the Poisson white noises, and the time delay in the interaction between the prey and predator is described approximately by a time delay parameter. The original ecosystem is first modeled as the Stratonovich SDE and then transferred to Itô SDE. Under some reasonable hypothesises, the

Conclusions
This paper investigates the stochastic behaviors of a time-delayed predator-prey ecosystem under pulse-type stochastic environment fluctuations. In this model, the stochastic fluctuations are characterized by the Poisson white noises, and the time delay in the interaction between the prey and predator is described approximately by a time delay parameter. The original ecosystem is first modeled as the Stratonovich SDE and then transferred to Itô SDE. Under some reasonable hypothesises, the stochastic averaging method is applied to simplify the original system. Finally, the approximate stationary PDFs of the prey and predator populations are obtained by using the perturbation method. To verify the accuracy of the proposed method, the results from Monte Carlo simulation are also calculated.
Based on the approximate stationary PDFs of prey and predator populations, the influences of the prey self-competition parameter, the time-delay parameter and the Poisson white noise excitation on the system are investigated in detail. The results show that, increasing time delay parameter as well as the mean arrival rate and the variance of the impulse of the Poisson white noise will enhance the fluctuations of the prey and predator population and make the ecosystem less stable. While the larger value of self-competition parameter will reduce the fluctuations of the system and increase the stability of the ecosystem.
In our present paper, only the effects of some system parameters are studied. Similar approaches may be adopted to study the influences of other system parameters. In the present investigation, the amplitudes of the Poisson white noises are assumed to be Gaussian distribution. They can be other types of distributions depending on the realistic situation.