Analysis of Stochastic Delay Predator-Prey System with Impulsive Toxicant Input in Polluted Environments

Abstract

A stochastic delay predator-prey model in a polluted environment with impulsive toxicant input is proposed and studied. The thresholds between stability in time average and extinction of each population are obtained. Some recent results are extended and improved greatly. Several simulation figures are introduced to support the conclusions.

1. Introduction

Environmental pollution by industries, agriculture, and other human activities is one of the most important socio-ecological problems in the world today. Due to toxins in the environment, lots of species have gone extinct, and many are on the verge of extinction. Thus, controlling the environmental pollution and the conservation of biodiversity are the major focus areas of all the countries around the world. This motivates scholars to study the effects of toxins on populations and to find out a theoretical persistence-extinction threshold.

Recently, a lot of population models in a polluted environment have been proposed and investigated; here, we may mention, among many others, [1–23]. Particularly, Yang et al. [15] pointed out that in many cases toxicants should be emitted in regular pulses, for example, the use of pesticides and the pollution by heavy metals (see, e.g., [24]). Thus, they proposed the following two-species Lotka-Volterra predator-prey system in a polluted environment with impulsive toxicant input: 𝑑𝑥₁(𝑡)𝑑𝑡=𝑥₁(𝑡)𝑟₁₀−𝑟₁₁𝐶₁₀(𝑡)−𝑎₁₁𝑥₁(𝑡)−𝑎₁₂𝑥₂(𝑡),𝑑𝑥₂(𝑡)𝑑𝑡=𝑥₂(𝑡)−𝑟₂₀−𝑟₂₁𝐶₂₀(𝑡)+𝑎₂₁𝑥₁(𝑡)−𝑎₂₂𝑥₂(𝑡),𝑑𝐶₁₀(𝑡)𝑑𝑡=𝑘₁𝐶_𝑒(𝑡)−𝑔₁+𝑚₁𝐶₁₀(𝑡),𝑑𝐶₂₀(𝑡)𝑑𝑡=𝑘₂𝐶_𝑒(𝑡)−𝑔₂+𝑚₂𝐶₂₀(𝑡),𝑑𝐶_𝑒(𝑡)𝑑𝑡=−ℎ𝐶_𝑒(𝑡),𝑡≠𝑛𝛾,𝑛∈𝑍⁺,Δ𝑥_𝑖(𝑡)=0,Δ𝐶_𝑖0(𝑡)=0,Δ𝐶_𝑒(𝑡)=𝑏,𝑡=𝑛𝛾,𝑛∈𝑍⁺,𝑖=1,2, (1) where all the parameters are positive constants andΔ𝑓(𝑡)=𝑓(𝑡⁺)−𝑓(𝑡), 𝑍⁺={1,2,…};𝑥₁(𝑡) and 𝑥₂(𝑡): the size of prey population and the predator population, respectively; 𝑟_𝑖0: the intrinsic growth rate of the 𝑖th population without toxicant; 𝑟_𝑖1: the 𝑖th population response to the pollutant present in the organism; 𝐶_𝑖0(𝑡): the concentration of toxicant in the 𝑖th organism; 𝐶_𝑒(𝑡): the concentration of toxicant in the environment; 𝑘𝐶_𝑒(𝑡): the organism’s net uptake of toxicant from the environment; 𝑔𝐶_𝑖0(𝑡)+𝑚𝐶_𝑖0(𝑡): the egestion and depuration rates of the toxicant in the 𝑖th organism; ℎ𝐶_𝑒(𝑡): the toxicant loss from the environment itself by volatilization and so on; 𝛾: the period of the impulsive effect about the exogenous input of toxicant; 𝑏: the toxicant input amount at every time.

Yang et al. [15] showed that in the following Lemma holds.

Lemma 1. For system (1), define Δ₂=𝑟₁₀𝑎₂₁−𝑟₂₀𝑎₁₁,Δ₂𝑎∶=²¹𝑟₁₁𝐾₁𝛾+𝑎₁₁𝑟₂₁𝐾₂𝛾,𝐾_𝑖=𝑘_𝑖𝑏ℎ𝑔_𝑖+𝑚_𝑖. (2)(a)If 𝑟₁₀<𝑟₁₁𝐾₁/𝛾, then lim_𝑡→+∞𝑥_𝑖(𝑡)=0, 𝑖=1,2.(b)If 𝑟₁₀>𝑟₁₁𝐾₁/𝛾 and Δ₂<Δ₂, then limsup_𝑡→+∞𝑡⁻¹∫^𝑡₀𝑥₁(𝑠)𝑑𝑠>0 and 𝑥₂(𝑡) goes to extinction.(c)If Δ₂>Δ₂, then limsup_𝑡→+∞𝑡⁻¹∫^𝑡₀𝑥_𝑖(𝑠)𝑑𝑠>0,𝑖=1,2.

Some interesting and important problems arise naturally.  (Q1) In the real world, the growth of species depends on various environmental factors, such as temperature, humidity and parasites and so forth. Therefore population models should be stochastic rather than deterministic (May [25]). Thus, what happens if model (1) is subject to stochastic noises? (Q2) In addition, time delays occur in almost every situation. Kuang [26] has pointed out that ignoring time delays means ignoring reality. Therefore, what happens if model (1) takes time delays into account? (Q3) Can we improve the results given in Lemma 1?

The aim of this paper is to study the above problems. Suppose that stochastic noises mainly affect the growth rates, with 𝑟_𝑖0→𝑟_𝑖0+𝛼_𝑖𝐵̇_𝑖(𝑡) (see, e.g., [27–39]), where 𝐵̇_𝑖(𝑡) is a white noise and 𝛼²_𝑖 is the intensity of the noise. Moreover, taking time delays into account, we obtain the following model: 𝑑𝑥₁(𝑡)=𝑥₁(𝑡)𝑟₁₀−𝑟₁₁𝐶₁₀(𝑡)−𝑎₁₁𝑥₁(𝑡)−𝑎₁₂𝑥₂𝑡−𝜏₁𝑑𝑡+𝛼₁𝑥₁(𝑡)𝑑𝐵₁(𝑡),𝑑𝑥₂(𝑡)=𝑥₂(𝑡)−𝑟₂₀−𝑟₂₁𝐶₂₀(𝑡)−𝑎₂₂𝑥₂(𝑡)+𝑎₂₁𝑥₁𝑡−𝜏₂𝑑𝑡+𝛼₂𝑥₂(𝑡)𝑑𝐵₂(𝑡),𝑑𝐶₁₀(𝑡)𝑑𝑡=𝑘₁𝐶_𝑒(𝑡)−𝑔₁+𝑚₁𝐶₁₀(𝑡),𝑑𝐶₂₀(𝑡)𝑑𝑡=𝑘₂𝐶_𝑒(𝑡)−𝑔₂+𝑚₂𝐶₂₀(𝑡),𝑑𝐶_𝑒(𝑡)𝑑𝑡=−ℎ𝐶_𝑒(𝑡),𝑡≠𝑛𝛾,𝑛∈𝑍⁺,Δ𝑥_𝑖(𝑡)=0,Δ𝐶_𝑖0(𝑡)=0,Δ𝐶_𝑒(𝑡)=𝑏,𝑡=𝑛𝛾,𝑛∈𝑍⁺,𝑖=1,2, (3) with initial condition 𝑥_𝑖(𝑡)=𝜙_𝑖(𝑡)>0,𝑡∈[−𝜏,0];𝜙_𝑖(0)>0,𝑖=1,2, (4) where 𝜏_𝑖≥0, 𝜏=max{𝜏₁,𝜏₂}, 𝜙_𝑖(𝑡) is continuous on [−𝜏,0]. Our main result is the following theorem.

Theorem 2. For system (3), define 𝜃₁=𝑟₁₀−0.5𝛼²₁,𝜃₂=𝑟₂₀+0.5𝛼²₂,Δ=𝑎₁₁𝑎₂₂+𝑎₁₂𝑎₂₁,Δ₁∶=𝜃₁𝑎₂₂+𝜃₂𝑎₁₂,Δ₂∶=𝜃₁𝑎₂₁−𝜃₂𝑎₁₁,Δ₁𝑎∶=²²𝑟₁₁𝐾₁𝛾−𝑎₁₂𝑟₂₁𝐾₂𝛾. (5)(i)If 𝜃₁<𝑟₁₁𝐾₁/𝛾, then both 𝑥₁ and 𝑥₂ go to extinction almost surely (a.s.); that is, lim_𝑡→+∞𝑥_𝑖(𝑡)=0𝑎.𝑠.,𝑖=1,2.(ii)If 𝜃₁>𝑟₁₁𝐾₁/𝛾 and Δ₂<Δ₂, then 𝑥₂(𝑡) goes to extinction and 𝑥₁ is stable in time average a.s.; that is, lim_𝑡→+∞𝑡^−1𝑡₀𝑥₁(𝑠)𝜃𝑑𝑠=¹−𝑟₁₁𝐾₁/𝛾𝑎₁₁>0,𝑎.𝑠. (6)(iii)If Δ₂>Δ₂, then both 𝑥₁ and 𝑥₂ are stable in time average a.s. lim_𝑡→+∞𝑡^−1𝑡₀𝑥₁(𝑠)Δ𝑑𝑠=¹−Δ₁Δ>0,𝑎.𝑠.lim_𝑡→+∞𝑡^−1𝑡₀𝑥₂(𝑠)Δ𝑑𝑠=²−Δ₂Δ>0,𝑎.𝑠. (7)

Remark 3. By comparing Lemma 1 with our Theorem 2, we can see that on the one hand, if 𝛼₁=𝛼₂=0 and 𝜏₁=𝜏₂=0, then 𝜃_𝑖=𝑟_𝑖0, Δ_𝑖=Δ_𝑖, 𝑖=1,2, and our stochastic delay system (3) becomes model (1); on the other hand, our results in Theorem 2 improve that in Lemma 1. Lemma 1 shows that the superior limit is positive, while Theorem 2 reveals that the limit exists and gives the explicit form of the limit. The contribution of this paper is therefore clear.

2. Proof

For the sake of simplicity, we introduce some notations: 𝑅²₊=𝑎=𝑎₁,𝑎₂∈𝑅²∣𝑎_𝑖,𝑓(𝑡>0,𝑖=1,2)=𝑡^−1𝑡₀𝑓(𝑠)𝑓(𝑡𝑑𝑠;)^∗=limsup_𝑡→+∞𝑡^−1𝑡₀𝑓(𝑠)𝑓(𝑡𝑑𝑠,)_∗=liminf_𝑡→+∞𝑡^−1𝑡₀𝑓(𝑠)𝑑𝑠. (8)

Lemma 4. For any given initial value 𝜙(𝑡)=(𝜙₁(𝑡),𝜙₂(𝑡))∈𝐶([−𝜏,0],𝑅²₊), there is a unique global positive solution 𝑥(𝑡)=(𝑥₁(𝑡),𝑥₂(𝑡))^𝑇 to the first two equations of system (3) a.s.

Proof. The proof is similar to Hung [29] by defining 𝑉(𝑥)=𝑥₁−1−ln𝑥₁+𝑎₁₁𝑎₂₂𝑎²₂₁𝑥₂−1−ln𝑥₂+0.5𝑎₁₁^𝑡_𝑡−𝜏2𝑥²₁(𝑠)𝑑𝑠 (9) and hence is omitted.

To begin with, let us consider the following subsystem of (3): 𝑑𝐶₁₀(𝑡)𝑑𝑡=𝑘₁𝐶_𝑒(𝑡)−𝑔₁+𝑚₁𝐶₁₀(𝑡),𝑑𝐶₂₀(𝑡)𝑑𝑡=𝑘₂𝐶_𝑒(𝑡)−𝑔₂+𝑚₂𝐶₂₀(𝑡),𝑑𝐶_𝑒(𝑡)𝑑𝑡=−ℎ𝐶_𝑒(𝑡),𝑡≠𝑛𝛾,𝑛∈𝑍⁺,Δ𝐶₁₀(𝑡)=0,Δ𝐶₂₀(𝑡)=0,Δ𝐶_𝑒(𝑡)=𝑏,𝑡=𝑛𝛾,𝑛∈𝑍⁺,0≤𝐶₀(0)≤1,0≤𝐶_𝑒(0)≤1. (10)

Lemma 5 (see [13, 15]). System (10) has a unique positive 𝛾-periodic solution (𝐶₁₀(𝑡),𝐶₂₀(𝑡),𝐶_𝑒(𝑡))^𝑇, and for each solution (𝐶₁₀(𝑡),𝐶₂₀(𝑡),𝐶_𝑒(𝑡))^𝑇 of (10), 𝐶₁₀(𝑡)→𝐶₁₀(𝑡), 𝐶₂₀(𝑡)→𝐶₂₀(𝑡), and 𝐶_𝑒(𝑡)→𝐶_𝑒(𝑡) as 𝑡→∞. Moreover, 𝐶_𝑖0(𝑡)>𝐶_𝑖0(𝑡) and 𝐶_𝑒(𝑡)>𝐶_𝑒(𝑡) for all 𝑡≥0 if 𝐶_𝑖0(0)>𝐶_𝑖0(0) and 𝐶_𝑒(0)>𝐶_𝑒(0), 𝑖=1,2, where 𝐶₁₀(𝑡)=𝐶₁₀(0)𝑒^{−(𝑔1+𝑚1)(𝑡−𝑛𝛾)}+𝑘₁𝑏𝑒^{−(𝑔1+𝑚1)(𝑡−𝑛𝛾)}−𝑒^{−ℎ(𝑡−𝑛𝛾)}ℎ−𝑔₁−𝑚₁1−𝑒^−ℎ𝛾,𝐶₂₀(𝑡)=𝐶₂₀(0)𝑒^{−(𝑔2+𝑚2)(𝑡−𝑛𝛾)}+𝑘₂𝑏𝑒^{−(𝑔2+𝑚2)(𝑡−𝑛𝛾)}−𝑒^{−ℎ(𝑡−𝑛𝛾)}ℎ−𝑔₂−𝑚₂1−𝑒^−ℎ𝛾,𝐶_𝑒(𝑡)=𝑏𝑒^{−ℎ(𝑡−𝑛𝛾)}1−𝑒^−ℎ𝛾,𝐶₁₀(0)=𝑘₁𝑏𝑒^{−(𝑔1+𝑚1)𝛾}−𝑒^−ℎ𝛾ℎ−𝑔₁−𝑚₁1−𝑒^{−(𝑔1+𝑚1)𝛾}1−𝑒^−ℎ𝛾,𝐶₂₀(0)=𝑘₂𝑏𝑒^{−(𝑔2+𝑚2)𝛾}−𝑒^−ℎ𝛾ℎ−𝑔₂−𝑚₂1−𝑒^{−(𝑔2+𝑚2)𝛾}1−𝑒^−ℎ𝛾,𝐶_𝑒(0)=𝑏1−𝑒^−ℎ𝛾 (11) for 𝑡∈(𝑛𝛾,(𝑛+1)𝛾] and 𝑛∈𝑍⁺. In addition, lim_𝑡→+∞𝑡^−1𝑡₀𝐶_𝑖0(𝑠)𝑘𝑑𝑠=^𝑖𝑏ℎ𝑔_𝑖+𝑚_𝑖𝛾=𝐾_𝑖𝛾,𝑖=1,2. (12)

Lemma 6 (see [34]). Suppose that 𝑥(𝑡)∈𝐶[Ω×[0,+∞),𝑅₊].  (I) If there exist 𝜎 and positive constants 𝜎₀,𝑇 such that ln𝑥(𝑡)≤𝜎𝑡−𝜎₀^𝑡₀𝑥(𝑠)𝑑𝑠+^𝑛_𝑖=1𝛽_𝑖𝐵_𝑖(𝑡) (13) for 𝑡≥𝑇, where 𝐵_𝑖(𝑡) are independent standard Brownian motions and 𝛽_𝑖 are constants, 1≤𝑖≤𝑛, then one has the following: if 𝜎≥0, then ⟨𝑥⟩^∗≤𝜎/𝜎₀ a.s.; if 𝜎<0, then lim_𝑡→+∞𝑥(𝑡)=0𝑎.𝑠.  (II) If there exist positive constants 𝜎₀,𝑇 and 𝜎 such that ln𝑥(𝑡)≥𝜎𝑡−𝜎₀^𝑡₀𝑥(𝑠)𝑑𝑠+^𝑛_𝑖=1𝛽_𝑖𝐵_𝑖(𝑡) (14) for 𝑡≥𝑇, then ⟨𝑥⟩_∗≥𝜎/𝜎₀ a.s.

Now, let us consider the following auxiliary system: 𝑑𝑦₁(𝑡)=𝑦₁(𝑡)𝑟₁₀−𝑟₁₁𝐶₁₀(𝑡)−𝑎₁₁𝑦₁(𝑡)𝑑𝑡+𝛼₁𝑦₁(𝑡)𝑑𝐵₁(𝑡),𝑑𝑦₂(𝑡)=𝑦₂(𝑡)−𝑟₂₀−𝑟₂₁𝐶₂₀(𝑡)+𝑎₂₁𝑦₁𝑡−𝜏₂−𝑎₂₂𝑦₂(𝑡)𝑑𝑡+𝛼₂𝑦₂(𝑡)𝑑𝐵₂(𝑡), (15) with initial value 𝜙(𝑡)∈𝐶([−𝜏,0],𝑅²₊).

Lemma 7. If 𝜃₁=𝑟₁₀−0.5𝛼²₁>𝑟₁₁𝐾₁/𝛾, then the solution 𝑦(𝑡) of system (15) obeys lim_𝑡→+∞𝑦₁(𝑡=𝜃)¹−𝑟₁₁𝐾₁/𝛾𝑎₁₁;lim_𝑡→+∞𝑦₂(𝑡)=0𝑎.𝑠.,𝑖𝑓Δ₂<Δ₂;lim_𝑡→+∞𝑦₂(𝑡=Δ)²−Δ₂𝑎₁₁𝑎₂₂𝑎.𝑠.,𝑖𝑓Δ₂>Δ₂. (16)

Proof. By Lemma 5, lim_𝑡→+∞𝑡^−1𝑡₀𝐶_𝑖0(𝑠)𝑑𝑠=lim_𝑡→+∞𝑡^−1𝑡₀𝐶_𝑖0(𝑠)=𝐾𝑑𝑠_𝑖𝛾,𝑖=1,2. (17) Then, for all𝜀>0, there exists 𝑇>0 such that 𝐾_𝑖𝛾−𝜀≤𝐶_𝑖0(𝑡≤𝐾)^𝑖𝛾+𝜀,𝑡>𝑇,𝑖=1,2. (18) An application of Itô’s formula to (15) yields ln𝑦₁(𝑡)−ln𝑦₁(0)=𝜃₁𝑡−𝑟₁₁^𝑡₀𝐶₁₀(𝑠)𝑑𝑠−𝑎₁₁^𝑡₀𝑦₁(𝑠)𝑑𝑠+𝛼₁𝐵₁(𝑡),ln𝑦₂(𝑡)−ln𝑦₂(0)=−𝜃₂𝑡−𝑟₂₁^𝑡₀𝐶₂₀(𝑠)𝑑𝑠+𝑎₂₁^𝑡₀𝑦₁𝑠−𝜏₂𝑑𝑠−𝑎₂₂^𝑡₀𝑦₂(𝑠)𝑑𝑠+𝛼₂𝐵₂(𝑡)=−𝜃₂𝑡−𝑟₂₁^𝑡₀𝐶₂₀(𝑠)𝑑𝑠+𝑎₂₁^𝑡₀𝑦₁(𝑠)𝑑𝑠−𝑎₂₁^𝑡_𝑡−𝜏2𝑦₁(𝑠)𝑑𝑠−⁰_−𝜏2𝑦₁(𝑠)𝑑𝑠−𝑎₂₂^𝑡₀𝑦₂(𝑠)𝑑𝑠+𝛼₂𝐵₂(𝑡). (19) That is to say, we have shown that 𝑡⁻¹𝑦ln¹(𝑡)𝑦₁(0)=𝜃₁−𝑟₁₁𝐶₁₀(𝑡)−𝑎₁₁𝑦₁(𝑡)+𝑡⁻¹𝛼₁𝐵₁(𝑡),𝑡⁻¹𝑦ln²(𝑡)𝑦₂(0)+𝑡⁻¹𝑎₂₁^𝑡_𝑡−𝜏2𝑦₁(𝑠)𝑑𝑠−⁰_−𝜏2𝑦₁(𝑠)𝑑𝑠=−𝜃₂−𝑟₂₁𝐶₂₀(𝑡)+𝑎₂₁𝑦₁(𝑡)−𝑎₂₂𝑦₂(𝑡)+𝑡⁻¹𝛼₂𝐵₂(𝑡). (20)(21) When (18) is used in (20), we can see that for 𝑡>𝑇, 𝑡⁻¹𝑦ln¹(𝑡)𝑦₁(0)≤𝜃₁−𝑟₁₁𝐾₁𝛾+𝑟₁₁𝜀−𝑎₁₁𝑦₁(𝑡)+𝑡⁻¹𝛼₁𝐵₁(𝑡),𝑡⁻¹𝑦ln¹(𝑡)𝑦₁(0)≥𝜃₁−𝑟₁₁𝐾₁𝛾−𝑟₁₁𝜀−𝑎₁₁𝑦₁(𝑡)+𝑡⁻¹𝛼₁𝐵₁(𝑡). (22)(23) Let 𝜀 be sufficiently small such that 𝜃₁−𝑟₁₁𝐾₁/𝛾−𝑟₁₁𝜀>0. Making use of (I) and (II) in Lemma 6 to (22) and (23), respectively, we have 𝑦₁(𝑡)^∗≤𝜃₁−𝑟₁₁𝐾₁/𝛾+𝑟₁₁𝜀𝑎₁₁,𝑦₁(𝑡)_∗≥𝜃₁−𝑟₁₁𝐾₁/𝛾−𝑟₁₁𝜀𝑎₁₁. (24) It then follows from the arbitrariness of 𝜀 that lim_𝑡→+∞𝑦₁(𝑡=𝜃)¹−𝑟₁₁𝐾₁/𝛾𝑎₁₁. (25) Substituting (17) and (25) into (20) and noting that lim_𝑡→+∞𝑡⁻¹𝐵₁(𝑡)=0, one can derive that lim_𝑡→+∞𝑡⁻¹ln𝑦₁(𝑡)=0,a.s. (26) Employing (20) and (21) in the expression 𝑎₂₁ln(𝑦₁(𝑡)/𝑦₁(0))+𝑎₁₁ln(𝑦₂(𝑡)/𝑦₂(0)) yields 𝑎₁₁𝑡⁻¹𝑦ln²(𝑡)𝑦₂(0)+𝑎₂₁𝑡⁻¹𝑦ln¹(𝑡)𝑦₁(0)=Δ₂−𝑟₁₁𝑎₂₁𝐶₁₀(𝑡)−𝑟₂₁𝑎₁₁𝐶₂₀(𝑡)−𝑎₁₁𝑎₂₂𝑦₂(𝑡)−𝑡⁻¹𝑎₁₁𝑎₂₁^𝑡_𝑡−𝜏2𝑦₁(𝑠)𝑑𝑠−⁰_−𝜏2𝑦₁(𝑠)𝑑𝑠+𝑡⁻¹𝑎₂₁𝛼₁𝐵₁(𝑡)+𝑎₁₁𝛼₂𝐵₂(𝑡). (27) In view of (25), we get lim_𝑡→+∞𝑡^−1𝑡_𝑡−𝜏2𝑦₁(𝑠)𝑑𝑠=lim_𝑡→+∞𝑡^−1𝑡₀𝑦₁(𝑠)𝑑𝑠−^𝑡−𝜏2₀𝑦₁(𝑠)𝑑𝑠=0,a.s. (28) By (17), (26), (27), and (28), for all𝜀>0, there exists 𝑇>0 such that, for 𝑡≥𝑇, 𝑎₁₁𝑡⁻¹𝑦ln²(𝑡)𝑦₂(0)≤Δ₂−Δ₂+𝜀−𝑎₁₁𝑎₂₂𝑦₂(𝑡)+𝑡⁻¹𝑎₂₁𝜎₁𝐵₁(𝑡)+𝑎₁₁𝜎₂𝐵₂(𝑡,𝑎₁₁)𝑡⁻¹𝑦ln²(𝑡)𝑦₂(0)≥Δ₂−Δ₂−𝜀−𝑎₁₁𝑎₂₂𝑦₂(𝑡)+𝑡⁻¹𝑎₂₁𝜎₁𝐵₁(𝑡)+𝑎₁₁𝜎₂𝐵₂(𝑡). (29)(30) If Δ₂<Δ₂, then we can choose 𝜀 sufficiently small such that Δ₂−Δ₂+𝜀<0. Then, by (29) and (I) in Lemma 6, we obtain lim_𝑡→+∞𝑦₂(𝑡)=0a.s. If Δ₂>Δ₂, then we can choose 𝜀 sufficiently small such that Δ₂−Δ₂−𝜀>0. An application of (I) and (II) in Lemma 6 to (29) and (30), respectively, makes one observe that Δ₂−Δ₂−𝜀𝑎₁₁𝑎₂₂≤𝑦₂(𝑡)_∗≤𝑦₂(𝑡)^∗≤Δ₂−Δ₂+𝜀𝑎₁₁𝑎₂₂,a.s. (31) Therefore, using the arbitrariness of 𝜀 results in lim_𝑡→+∞𝑦₂(𝑡=Δ)²−Δ₂𝑎₁₁𝑎₂₂a.s. (32) This completes the proof.

We are now in the position to prove our main results.

Proof of Theorem 2. Applying Itô’s formula to (3) leads to ln𝑥₁(𝑡)−ln𝑥₁(0)=𝜃₁𝑡−𝑟₁₁^𝑡₀𝐶₁₀(𝑠)𝑑𝑠−𝑎₁₁^𝑡₀𝑥₁(𝑠)𝑑𝑠−𝑎₁₂^𝑡₀𝑥₂𝑠−𝜏₁𝑑𝑠+𝛼₁𝐵₁(𝑡)=𝜃₁𝑡−𝑟₁₁^𝑡₀𝐶₁₀(𝑠)𝑑𝑠−𝑎₁₂^𝑡₀𝑥₂(𝑠)𝑑𝑠+𝑎₁₂^𝑡_𝑡−𝜏1𝑥₂(𝑠)𝑑𝑠−⁰_−𝜏1𝑥₂(𝑠)𝑑𝑠−𝑎₁₁^𝑡₀𝑥₁(𝑠)𝑑𝑠+𝛼₁𝐵₁(𝑡).ln𝑥₂(𝑡)−ln𝑥₂(0)=−𝜃₂𝑡−𝑟₂₁^𝑡₀𝐶₂₀(𝑠)𝑑𝑠+𝑎₂₁^𝑡₀𝑥₁(𝑠)𝑑𝑠−𝑎₂₁^𝑡_𝑡−𝜏2𝑥₁(𝑠)𝑑𝑠−⁰_−𝜏2𝑥₁(𝑠)𝑑𝑠−𝑎₂₂^𝑡₀𝑥₂(𝑠)𝑑𝑠+𝛼₂𝐵₂(𝑡). (33)(34)
(i) It follows from (17) and (33) that 𝑡⁻¹ln𝑥₁(𝑡)−𝑡⁻¹ln𝑥₁(0)≤𝜃₁−𝑟₁₁𝐶₁₀(𝑡)−𝑎₁₁𝑥₁(𝑡+𝛼)¹𝐵₁(𝑡)𝑡≤𝜃₁−𝑟₁₁𝐾₁𝜆+𝜀−𝑎₁₁𝑥₁(𝑡+𝛼)¹𝐵₁(𝑡)𝑡 (35) for sufficiently large 𝑡. Since 𝜃₁−𝑟₁₁𝐾₁/𝜆<0, then we can choose 𝜀 sufficiently small such that 𝜃₁−𝑟₁₁𝐾₁/𝜆+𝜀<0. Then, by (I) in Lemma 6, lim_𝑡→+∞𝑥₁(𝑡)=0,a.s. (36) When (36) is used in (34), one can see that 𝑡⁻¹ln𝑥₂(𝑡)−ln𝑥₂(0)≤−𝜃₂+𝜀−𝑎₂₂𝑥₂(𝑡+𝛼)²𝐵₂(𝑡)𝑡 (37) for sufficiently large 𝑡, where 𝜀>0 obeys −𝜃₂+𝜀<0. In view of Lemma 6 again, lim_𝑡→+∞𝑥₂(𝑡)=0,a.s.
(ii) By the stochastic comparison theorem [40], one can observe that 𝑥₁(𝑡)≤𝑦₁(𝑡),𝑥₂(𝑡)≤𝑦₂(𝑡). (38) Note that 𝜃₁>𝑟₁₁𝐾₁/𝛾 and Δ₂<Δ₂; it then follows from Lemma 7 that lim_𝑡→+∞𝑦₂(𝑡)=0,a.s. Making use of (38) gives lim_𝑡→+∞𝑥₂(𝑡)=0,a.s. Thus, for all𝜀>0, there exists 𝑇>0 such that, for 𝑡≥𝑇, 𝜀2≤𝑎₁₂𝑥₂(𝑡)≤𝜀2. (39) Substituting the above inequalities into (33) and then using (18), we obtain 𝑡⁻¹ln𝑥₁(𝑡)≤𝑡⁻¹ln𝑥₁(0)+𝜃₁−𝑟₁₁𝐶₁₀(𝑡)−𝑎₁₁𝑥₁(𝑡+𝜀)2+𝛼₁𝐵₁(𝑡)𝑡≤𝜃₁−𝑟₁₁𝐾₁𝛾+2𝜀−𝑎₁₁𝑥₁(𝑡+𝛼)¹𝐵₁(𝑡)𝑡,𝑡⁻¹ln𝑥₁(𝑡)≥𝑡⁻¹ln𝑥₁(0)+𝜃₁−𝑟₁₁𝐶₁₀(𝑡)−𝑎₁₁𝑥₁(𝑡−𝜀)2+𝛼₁𝐵₁(𝑡)𝑡≥𝜃₁−𝑟₁₁𝐾₁𝛾−2𝜀−𝑎₁₁𝑥₁(𝑡+𝛼)¹𝐵₁(𝑡)𝑡. (40)(41) Let 𝜀 be sufficiently small such that 𝜃₁−𝑟₁₁𝐾₁/𝛾−𝜀>0, and then, applying (I) and (II) in Lemma 6 to (40) and (41), respectively, one can see that 𝜃₁−𝑟₁₁𝐾₁/𝛾−2𝜀𝑎₁₁≤𝑥₁(𝑡)_∗≤𝑥₁(𝑡)^∗≤𝜃₁−𝑟₁₁𝐾₁/𝛾+2𝜀𝑎₁₁a.s. (42) An application of the arbitrariness of 𝜀 gives lim_𝑡→+∞𝑥₁(𝑡=𝜃)¹−𝑟₁₁𝐾₁/𝛾𝑎₁₁,a.s. (43)
(iii) Clearly, Δ₂>Δ₂ implies 𝜃₁>𝑟₁₁𝐾₁/𝛾, and then, by Lemma 7, lim_𝑡→+∞𝑦₂(𝑡=Δ)²−Δ₂𝑎₁₁𝑎₂₂. (44) Thus, similar to the proof of (28), we get lim_𝑡→+∞𝑡^−1𝑡_𝑡−𝜏1𝑦₂(𝑠)𝑑𝑠=0a.s. (45) Therefore, by (26), (28), and (38), we can observe that limsup_𝑡→+∞𝑡⁻¹ln𝑥₁(𝑡)≤lim_𝑡→+∞𝑡⁻¹ln𝑦₁(𝑡)=0,lim_𝑡→+∞𝑡^−1𝑡_𝑡−𝜏2𝑥₁(𝑠)𝑑𝑠=0,lim_𝑡→+∞𝑡^−1𝑡_𝑡−𝜏1𝑥₂(𝑠)𝑑𝑠=0,a.s. (46)(47) Employing (33) and (34) in the expression 𝑎₂₁ln(𝑥₁(𝑡)/𝑥₁(0))+𝑎₁₁ln(𝑥₂(𝑡)/𝑥₂(0)) yields 𝑡⁻¹𝑎₂₁𝑥ln¹(𝑡)𝑥₁(0)+𝑡⁻¹𝑎₁₁𝑥ln²(𝑡)𝑥₂(0)=𝑎₁₂𝑎₂₁𝑡^−1𝑡_𝑡−𝜏1𝑥₂(𝑠)𝑑𝑠−⁰_−𝜏1𝑥₂(𝑠)𝑑𝑠−𝑎₁₁𝑎₂₁𝑡^−1𝑡_𝑡−𝜏2𝑥₁(𝑠)𝑑𝑠−⁰_−𝜏2𝑥₁(𝑠)+Δ𝑑𝑠₂−𝑎₂₁𝑟₁₁𝐶₁₀(𝑡)−𝑎₁₁𝑟₂₁𝐶₂₀(𝑡𝑥)−Δ₂(𝑡)+𝑡⁻¹𝑎₂₁𝛼₁𝐵₁(𝑡)+𝑡⁻¹𝑎₁₁𝛼₂𝐵₂(𝑡). (48) When (18), (46) and (47), are used in (48), one can obtain 𝑡⁻¹𝑎₁₁𝑥ln²(𝑡)𝑥₂(0)≥Δ₂−Δ₂−𝜀−Δ𝑥₂(𝑡)+𝑡⁻¹𝑎₂₁𝛼₁𝐵₁(𝑡)+𝑡⁻¹𝑎₁₁𝛼₂𝐵₂(𝑡) (49) for sufficiently large 𝑡, where 𝜀>0 obeys Δ₂−Δ₂−𝜀>0. It then follows from (II) in Lemma 6 that 𝑥₂(𝑡)_∗≥Δ₂−Δ₂−𝜀Δ. (50) By virtue of the arbitrariness of 𝜀, we can see that 𝑥₂(𝑡)_∗≥Δ₂−Δ₂Δ. (51) Consequently, for every 0<𝜀<𝑎₁₂(Δ₂−Δ₂)/Δ, there is 𝑇>0 such that 𝑎₁₂𝑥₂(𝑡)≥𝑎₁₂𝑥_2∗𝑎−𝜀≥¹²Δ₂−Δ₂Δ−𝜀,𝑡>𝑇. (52) Substituting the above inequality into (33) and then using (18) and (47), one can see that 𝑡⁻¹𝑥ln¹(𝑡)𝑥₁(0)≤𝜃₁−𝑎₁₂Δ₂−Δ₂Δ+3𝜀−𝑎₁₁𝑡^−1𝑡₀𝑥₁(𝑠)𝑑𝑠+𝑡⁻¹𝛼₁𝐵₁(𝑡)=𝑎₁₁Δ₁−Δ₁Δ+3𝜀−𝑎₁₁𝑡^−1𝑡₀𝑥₁(𝑠)𝑑𝑠+𝑡⁻¹𝛼₁𝐵₁(𝑡) (53) for sufficiently large 𝑡. Since Δ₁−Δ₁>0, and then, by Lemma 6 and the arbitrariness of 𝜀, one can observe that 𝑥₁(𝑡)^∗≤Δ₁−Δ₁Δ. (54) When this inequality, (18) and (47), are used in (34), we can see that 𝑡⁻¹𝑥ln²(𝑡)𝑥₂(0)≤−𝜃₂+𝑎₂₁Δ₁−Δ₁Δ+3𝜀−𝑎₂₂𝑡^−1𝑡₀𝑥₂(𝑠)𝑑𝑠+𝑡⁻¹𝛼₂𝐵₂(𝑡)=𝑎₂₂Δ₂−Δ₂Δ+3𝜀−𝑎₂₂𝑡^−1𝑡₀𝑥₂(𝑠)𝑑𝑠+𝑡⁻¹𝛼₂𝐵₂(𝑡) (55) for sufficiently large 𝑡. Then, it follows from Lemma 6 and the arbitrariness of 𝜀 that 𝑥₂(𝑡)^∗≤Δ₂−Δ₂Δ. (56) Substituting the above inequality and (18) into (33), we get 𝑡⁻¹𝑥ln¹(𝑡)𝑥₁(0)≥𝜃₁−𝑎₁₂Δ₂−Δ₂Δ−3𝜀−𝑎₁₁𝑡^−1𝑡₀𝑥₁(𝑠)𝑑𝑠+𝑡⁻¹𝛼₁𝐵₁(𝑡)=𝑎₁₁Δ₁−Δ₁Δ−3𝜀−𝑎₁₁𝑡^−1𝑡₀𝑥₁(𝑠)𝑑𝑠+𝑡⁻¹𝛼₁𝐵₁(𝑡) (57) for sufficiently large 𝑡. By (II) in Lemma 6 and the arbitrariness of 𝜀 again, we obtain 𝑥₁(𝑡)_∗≥Δ₁−Δ₁Δ. (58) Then, the required assertion follows from (51), (54), (56), and (58).

3. Numerical Simulations

Let us use the famous Milstein method (see, e.g., [41]) to illustrate the analytical results.

To begin with, we choose 𝑟₁₀=0.85, 𝑟₂₀=0.05, 𝑟₁₁=𝑟₂₁=1, 𝑎₁₁=0.4,   𝑎₁₂=0.4, 𝑎₂₁=0.3, 𝑎₂₂=0.3, 𝜏₁=3,𝜏₂=8, 𝛼²₂=0.1, 𝑘_𝑖=𝑔_𝑖=𝑚_𝑖=0.1, 𝑖=1,2, ℎ=0.5, 𝑏=0.6, and 𝛾=12. Then, 𝐾_𝑖=𝑘_𝑖𝑏ℎ𝑔_𝑖+𝑚_𝑖Δ=0.6,₂=𝑟₁₀𝑎₂₁−𝑟₂₀𝑎₁₁=0.235>Δ₂=𝑎₂₁𝑟₁₁𝐾₁𝛾+𝑎₁₁𝑟₂₁𝐾₂𝛾=0.035. (59) By (c) in Lemma 1, the solution of model (1) obeys limsup_𝑡→+∞𝑡^−1𝑡₀𝑥₁(𝑠)𝑑𝑠>0,limsup_𝑡→+∞𝑡^−1𝑡₀𝑥₂(𝑠)𝑑𝑠>0. (60) However, when the white noises are taken into account, the properties of the system may be changed greatly. In Figure 1, we let the coefficients be same with the above. The only difference between conditions of Figures 1(a), 1(b), and 1(c) is that the value of 𝛼²₁ is different. In Figure 1(a), we choose 𝛼²₁/2=0.82. Therefore, 𝜃₁=𝑟₁₀−𝛼²₁2𝑟=0.03<¹¹𝐾₁𝛾=0.05. (61) Then, by (i) in Theorem 2, both 𝑥₁ and 𝑥₂ are extinctive. Figure 1(a) confirms these. In Figure 1(b), we choose 𝛼²₁/2=0.65. That is to say 𝜃₁=0.2>𝑟₁₁𝐾₁/𝛾=0.05 and Δ₂=0.02<Δ₂=0.035. It then follows from (ii) in Theorem 2 that 𝑥₂ is extinctive and 𝑥₁ is stable in time average: lim_𝑡→+∞𝑥₁(𝑡=𝜃)¹−𝑟₁₁𝐾₁/𝛾𝑎₁₁=0.375. (62) See Figure 1(b). In Figure 1(c), we choose 𝛼²₁/2=0.2. Then, Δ₂=0.155>Δ₂=0.035. In view of (iii) in Theorem 2, we can obtain that both 𝑥₁ and 𝑥₂ are stable in time average: lim_𝑡→+∞𝑥₁(𝑡=Δ)¹−Δ₁Δ=0.240.24=1,lim_𝑡→+∞𝑥₂(𝑡=Δ)²−Δ₂Δ=0.120.24=0.5. (63) Figure 1(c) confirms these.

(a)

(b)

(c)

Figure 1 

Solutions of system (3) for 𝑟10=0.85, 𝑟20=0.05, 𝑟11=𝑟21=1, 𝑎11=0.4, 𝑎12=0.4, 𝑎21=0.3, 𝑎22=0.3, 𝜏1=3, 𝜏2=8, 𝛼22=0.1, 𝑘𝑖=𝑔𝑖=𝑚𝑖=0.1, 𝑖=1,2, ℎ=0.5, 𝑏=0.6, 𝛾=12, 𝑥1(0)=0.9, 𝑥2(0)=0.5, 𝐶0(0)=𝐶𝑒(0)=0.1, and step size Δ𝑡=0.001. (a) is with 𝛼21/2=0.82; (b) is with 𝛼21/2=0.65; (c) is with 𝛼21/2=0.2.

In Figure 2, we choose 𝑟₁₀=0.85, 𝑟₂₀=0.05, 𝑟₁₁=𝑟₂₁=1, 𝑎₁₁=0.4, 𝑎₁₂=0.4, 𝑎₂₁=0.3, 𝑎₂₂=0.3, 𝜏₁=3, 𝜏₂=8, 𝛼²₁=0.4, 𝛼²₂=0.1, 𝑘_𝑖=𝑔_𝑖=𝑚_𝑖=0.1, 𝑖=1,2, ℎ=0.5, and 𝑏=0.6. The only difference between conditions of Figures 1(c) and 2 is that the value of 𝛾 is different. In Figure 2, we choose 𝛾=0.8. Then, 𝜃₁=0.65<𝑟₁₁𝐾₁/𝛾=0.75. It follows from (i) in Theorem 2 that both 𝑥₁ and 𝑥₂ are extinctive. Figure 2 confirms these. By comparing Figure 1(c) with Figure 2, one can see that the impulsive period 𝛾 plays a key role in determining the stability in time average and the extinction of the species.

Figure 2 

Solutions of system (3) for 𝑟10=0.85, 𝑟20=0.05, 𝑟11=𝑟21=1, 𝑎11=0.4, 𝑎12=0.4, 𝑎21=0.3, 𝑎22=0.3, 𝜏1=3, 𝜏2=8, 𝛼21=0.4, 𝛼22=0.1, 𝑘𝑖=𝑔𝑖=𝑚𝑖=0.1, 𝑖=1,2, ℎ=0.5, 𝑏=0.6, 𝛾=0.8, 𝑥1(0)=0.9, 𝑥2(0)=0.5, 𝐶0(0)=𝐶𝑒(0)=0.1, and step size Δ𝑡=0.001.

4. Conclusions and Future Directions

This paper is concerned with stochastic delay predator-prey model in a polluted environment with impulsive toxicant input. For each species, the threshold between stability in time average and extinction is established. Some recent results are improved and extended. Our Theorem 2 reveals some interesting and important results. (A)Firstly, time delay is harmless for stability in time average and extinction of the stochastic system (3).(B)The white noise 𝛼₁𝑑𝐵₁(𝑡) and 𝛼₂𝑑𝐵₂(𝑡) can change the properties of the system greatly.(C)The impulsive period 𝛾 plays an important role in determining the stability in time average and the extinction of the species.

Some interesting questions deserve further investigations. One may consider some more realistic but more complex systems, for example, stochastic delay model with Markov switching (see, e.g., [30, 32, 39]). It is also interesting to investigate what happens if 𝑎_𝑖𝑗 is stochastic.

Acknowledgments

The author thanks the editor and reviewer for these valuable and important comments. This research is supported by NSFC of China (nos. 11301207, 11171081, 11301112 and 11171056), Natural Science Foundation of Jiangsu Province (No. BK20130411) and Natural Science Research Project of Ordinary Universities in Jiangsu Province (no. 13KJB110002).