Electronic Journal of Differential Equations, Vol. 2022 (2022), No. 32, pp. 1–28. ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu DYNAMICS OF STOCHASTIC LOTKA-VOLTERRA PREDATOR-PREY MODELS DRIVEN BY THREE INDEPENDENT BROWNIAN MOTIONS SHANGZHI LI, SHANGJIANG GUO Abstract. This article concerns the permanence and extinction of stochastic Lotka-Volterra predator-prey models perturbed by three independent white noises. We establish some criteria and present some numerical simulations that illustrate our theoretical results. It is shown that the presence of strong noise on either the intra-specific interaction rate or the inter-specific interaction rate may lead to complete different dynamical behaviors from the deterministic case. 1. Introduction Predator-prey interaction is one of the basic interspecies relations in nature and society [3, 6, 11, 26]. Because of environmental stochastic perturbations, more and more researchers are interest on the following stochastic predator-prey model driven by three independent Brownian motions W1(·), W2(·) and W3(·): dX1(t) = X1(t)(a1 − b1X1(t)− c1X2(t))dt + (µ1X 2 1 (t) + γ1X1(t))dW1(t) + ρ1X1(t)X2(t)dW3(t), dX2(t) = X2(t)(a2 − b2X2(t) + c2X1(t))dt + (µ2X 2 2 (t) + γ2X2(t))dW2(t) + ρ2X1(t)X2(t)dW3(t). (1.1) where X1(t) and X2(t) denote the quantities of prey and predator populations, respectively, a1 and a2 are intrinsic growth rates, positive constants b1 and b2 represent the intra-specific interaction rates, positive constants c1 and c2 represent the inter-specific interaction. More precisely, c1 is the death rate per encounter of prey due to predation, and c2/c1 is the efficiency of turning predated prey into predator. The dynamics of (1.1) without noise (i.e., µj = γj = ρj = 0 for j = 1, 2) is quite simple [1, 22]. However, the presence of noise makes (1.1) have more complicated dynamical behaviour. For example, in [2, 9, 30, 32]) there are a plenty of results on the stochastic predator-prey model dX1(t) = X1(t)(a1 − b1X1(t)− c1X2(t))dt+ γ1X1(t)dW (t), dX2(t) = X2(t)(−a2 − b2X2(t) + c2X1(t))dt+ γ2X2(t)dW (t), (1.2) 2020 Mathematics Subject Classification. 34C12, 60H10, 92D25. Key words and phrases. Predator-prey model; Lyapunov exponent; permanence; extinction. ©2022. This work is licensed under a CC BY 4.0 license. Submitted March 21, 2021. Published April 20, 2022. 1 2 S. LI, S. GUO EJDE-2022/32 where all the parameters are positive and the intrinsic growth rates are perturbed stochastically. Arnold, Horsthemke and Stucki [2] investigated the sample paths of equation (1.2). Rudniki et al [27, 28, 29] investigated the convergence of densities of the distributions of the solutions of (1.2). Mao, Sabanis and Renshaw [21] investigated the existence and uniqueness of the positive solution of the following model, where the intra-specific interaction rates are perturbed stochastically, dX1(t) = X1(t)(a1 − b1X1(t)− c1X2(t))dt+ µ1X 2 1 (t)dW (t), dX2(t) = X2(t)(−a2 − b2X2(t) + c2X1(t))dt+ µ2X 2 2 (t)dW (t). (1.3) Dang, Du, and Ton [5] proved that the densities of (1.3) either converge in L1 to an invariant density or converge weakly to a singular measure on the boundary as well. It is natural to ask what happens to dynamical behaviour of model (1.1) in which there are noises and perturbations on both inter-specific interaction rates intra- specific interaction rates. As we know, there are few literatures in this direction. For convenience, by assuming γ1 = γ2 = 0 we shall investigate the system dX1(t) = X1(t)(a1 − b1X1(t)− c1X2(t))dt + µ1X 2 1 (t)dW1(t) + ρ1X1(t)X2(t)dW3(t), dX2(t) = X2(t)(a2 − b2X2(t) + c2X1(t))dt + µ2X 2 2 (t)dW2(t) + ρ2X1(t)X2(t)dW3(t). (1.4) Stochastic perturbations on intrinsic growth rates a1 and a2 will be investigated in our another paper. Much effort has been devoted to the study of prey-predator systems with a1 > 0 > a2 (that is, the intrinsic growth rate of the prey population is positive while the intrinsic growth rate of the predator population is negative); See, for example, [5, 10, 16, 21]. In such prey-predator systems, the predator population relies on a single species for food. In a real world, most species feed on more than one species. Therefore, in this paper we shall not confine our attention only to the case where a1 > 0 > a2 and shall distinguish four cases (i.e., a1 < 0 and a2 < 0, a1 < 0 < a2, a1 > 0 > a2, a1 > 0 and a2 > 0) to investigate the dynamics of system (1.4). To make the ecological model more accurate, we consider three different white noises in one model, since both the intra-specific interaction rates and inter-specific interaction rates might be perturbed by environmental randomness. One of important concepts in stochastic population models is stochastic perma- nence, which indicates that the species will survive forever. Thus, another purpose is to describe the effect of the three white noises on the stochastic permanence of (1.4). In particular, we want to see whether and how large intensities of noises could lead to extinction even though the population persists in the associated de- terministic system, and also to see whether and how large intensities of noises could lead to the permanence of the population even though some population dies out in the associated deterministic system. Finding conditions ensuring the stochastic permanence is drawing plenty of at- tentions (see, for example, [10]). Some useful methods such as Lyapunov-type functions and ergodicity have been proposed [8, 17]. Thus, it is interesting to develop much sharper or more general criteria by formulating more general and better candidates for Lyapunov functions. Recently, Nguyen and Yin [23] studied the coexistence and exclusion of stochastic competitive Lotka-Volterra models by establishing a threshold in terms of dynamics on the boundary. Different from EJDE-2022/32 STOCHASTIC LOTKA-VOLTERRA PREDATOR-PREY MODELS 3 the Lyapunov method, in this paper we shall develop Nguyen and Yin’s approach in [23] and introduce two thresholds λ1 and λ2 to determine the permanence and extinction in stochastic prey-predator system (1.4). In particular, we know that Lotka-Volterra models laid a theoretical foundation for competition among species [13, 18, 19]. Here, we should point out that our approach is applicable for a variety of other stochastic Lotka-Volterra competition models. Let (Ω,F , {F}t≥0,P) be a complete filtered probability space with the filtration {F}t≥0. We suppose that µi 6= 0, i = 1, 2, so that the diffusion is non-degenerate. LetXx(t) = (X1,x(t), X2,x(t)) be the solution to (1.4) with initial value x = (x1, x2). Denote that R2,◦ + = {(x1, x2) : x1 > 0, x2 > 0}. Throughout this paper, suppose that µi > 0, ρi ≥ 0, i = 1, 2. For each x ∈ R2,◦ + , it is proved in [20] that the solution Xx(t) remains in R2,◦ + with probability 1. It is easy to see that the solution Xx(t) is a Markov process. This article is organized as follows. In Section 2, we investigate the dynamics on the boundary of the solutions and then derive some thresholds that are used to determine extinction and permanence. In particular, we state our main theoret- ical results and describe the main difference between the stochastic system (1.4) and its associated deterministic system. In Section 3, we perform some numerical simulations to illustrate our theoretical results. Section 4 is devoted to the proof of the extinction and weak convergence to a boundary distribution of species of system (1.4) with a1 < 0. Section 5 is devoted to the dynamics of system (1.4) with a1 > 0 > a2. In Section 5, we provide the proof of the weak convergence to a boundary distribution and coexistence of the species of (1.4) in the case where a1 > 0 and a2 > 0. 2. Main results To state our main result, we need to introduce an auxiliary process. For each fixed j = 1, 2, consider dψ(t) = ψ(t)(aj − bjψ(t))dt+ µjψ 2(t)dWj(t), (2.1) where W1(t) and W2(t) are defined as in (1.4). Let ψj,x be the solution to (2.1) starting at xj . Equation (2.1) with aj > 0 has a unique invariant probability measure π∗j in (0,∞) with density (see [5] for more details) f∗j (φ) = c∗j φ4 exp ( 2bj µ2 jφ − aj µ2 jφ 2 ) for φ > 0. Here c∗j is a normalizing constant. Moreover, for every x > 0 and p < 3, P { lim T→∞ 1 T ∫ T 0 ψpj,x(t)dt = Qjp , ∫ ∞ 0 φpf∗j (φ)dφ } = 1. (2.2) By direct calculation, we have aj = bjQj1 + µ2 j 2 Qj2 when aj > 0. (2.3) The following two quantities play an important role in our analysis: λ1 = a2 + c2Q11 − ρ2 2 2 Q12 when a1 > 0, λ2 = a1 − c1Q21 − ρ2 1 2 Q22 when a2 > 0. 4 S. LI, S. GUO EJDE-2022/32 We start with the case where the prey population has a negative intrinsic growth rate, and have the following result on the stochastic system (1.4) with a1 < 0. Theorem 2.1. (i) If a1 < 0 and a2 < 0 then every solution Xx(t) to (1.4) with initial value x ∈ R2,◦ + satisfies that almost surely limt→∞Xx(t) = (0, 0). (ii) If a1 < 0 < a2 then every solution Xx(t) = (X1,x(t), X2,x(t)) to (1.4) with initial value x ∈ R2,◦ + satisfies that almost surely limt→∞X1,x(t) = 0 and the distribution of X2,x(t) converges weakly to π∗2 . As for the deterministic system associated with (1.4), we see that both prey and predator population die out when a1 < 0 and a2 < 0, and that the prey population goes extinct and predator tends to a2/b2 when a1 < 0 < a2. Comparing with the deterministic system associated with (1.4), Theorem 2.1 implies that the three white noises affect only the the amplitude of the oscillation of sample paths instead of the long-time behavior of solutions to system (1.4) with a1 < 0. In particular, the assumption that a1 < 0 implies that almost surely limt→∞ ψ1,x(t) = 0, and hence the density of the unique invariant probability measure of (2.1) is exactly the Dirac delta function δ(·). In this case, λ1 is actually equal to a2. As we shall see later, the sign of λ1 determines whether the predator dies out or not. Theorem 2.2. (i) If a1 > 0 > a2 and λ1 < 0 then every solution Xx(t) = (X1,x(t), X2,x(t)) to (1.4) with initial value x ∈ R2,◦ + satisfies that almost surely limt→∞X2,x(t) = 0 and X1,x(t) converges weakly to π∗1 . (ii) If a1 > 0 > a2 and λ1 > 0 then for any initial point x ∈ R2,◦ + , system (1.4) is permanent, i.e., its solution Xx(t) has a unique invariant probability concentrated on R2,◦ + . It is easy to see that Theorem 2.2 unifies and improves the results obtained by Dang et al [5] who investigated the asymptotic behavior of system (1.3). As for the deterministic system associated with (1.4), we see that there exists a unique globally asymptotically stable (GAS) equilibrium (a1/b1, 0) when 0 < a1c2 < −a2b1, and that there exists a unique GAS equilibrium ((a1b2 − a2c1)/(b1b2 + c1c2), (a1c2 + a2b1)/(b1b2 + c1c2)) when a1c2 > −a2b1 > 0. In particular, if 0 < a1c2 < −a2b1, then it follows from (2.3) that λ1 < a2 − a2b1 a1 Q11 − ρ2 2 2 Q12 = a2µ 2 1 2a1 Q12 − ρ2 2 2 Q12 < 0. In view of Theorem 2.2(i), we obtain the following result. Corollary 2.3. Assume that 0 < a1c2 < −a2b1, then every solution Xx(t) = (X1,x(t), X2,x(t)) to (1.4) with initial value x ∈ R2,◦ + satisfies that almost surely limt→∞X2,x(t) = 0, X1,x(t) converges weakly to π∗1 . The above corollary implies that if the predator population of the deterministic system associated with (1.4) dies out eventually, then the presence of noise in (1.4) cannot change this extinction tendency but either accelerate or delay the extinction of the predator population. However, if there is a stable coexisting equilibrium in the deterministic system associated with (1.4), then Theorem 2.2 implies that the presence of immoderate noise may cause the predator population to die out as well. This is exactly the essential difference between the stochastic system (1.4) and its deterministic system. This interesting observation can also be obtained in the case where a1 > 0 and a2 > 0. Namely, we have the following result. EJDE-2022/32 STOCHASTIC LOTKA-VOLTERRA PREDATOR-PREY MODELS 5 Theorem 2.4. Assume that a1 > 0 and a2 > 0. (i) If λ1 > 0 and λ2 < 0 then every solution Xx(t) = (X1,x(t), X2,x(t)) to (1.4) with initial value x ∈ R2,◦ + satisfies that almost surely lim supt→∞X2,x(t) > 0, limt→∞X1,x(t) = 0 and X2,x(t) converges weakly to π∗2 . (ii) If λ1 < 0 and λ2 > 0 then every solution Xx(t) = (X1,x(t), X2,x(t)) to (1.4) with initial value x ∈ R2,◦ + satisfies that almost surely lim supt→∞X1,x(t) > 0, limt→∞X2,x(t) = 0 and X1,x(t) converges weakly to π∗1 . (iii) If λ1 > 0 and λ2 > 0 then for any initial point x ∈ R2,◦ + , then system (1.4) is permanent, i.e., its solution Xx(t) has a unique invariant probability concentrated on R2,◦ + . If a2 > 0 and ρ2 = 0 (i.e., there is no random influence on the inter-specific interaction term of prey population), then we have λ1 > 0, which together with Theorem 2.4(i)(iii) implies the permanence of X2,x(t) for all x ∈ R2,◦ + and the existence of an invariant probability measure. Namely, we have the following result. Corollary 2.5. Assume that a1 > 0, a2 > 0, and ρ2 = 0. (i) If λ2 < 0 then every solution Xx(t) = (X1,x(t), X2,x(t)) to (1.4) with ini- tial value x ∈ R2,◦ + satisfies that almost surely lim supt→∞X2,x(t) > 0, limt→∞X1,x(t) = 0 and X2,x(t) converges weakly to π∗2 . (ii) If λ2 > 0 then for any initial point x ∈ R2,◦ + , then system (1.4) is permanent, i.e., its solution Xx(t) has a unique invariant probability concentrated on R2,◦ + . The dynamics of the deterministic system associated with (1.4) with a1 > 0 and a2 > 0 can be classified in terms of the sign of a1b2 − a2c1. In fact, the associated system has a unique GAS equilibrium (0, a2/b2) when 0 < a1b2 < a2c1, and has a unique GAS equilibrium ((a1b2−a2c1)/(b1b2+c1c2), (a1c2+a2b1)/(b1b2+ c1c2)) when a1b2 > a2c1 > 0. Even though the associated deterministic system is permanent, as we shall see, the presence of strong noises on the intra-specific interaction rate of the prey and predator may result in the extinction of the prey or predator population. On the other hand, even though the prey population of the associated deterministic system goes to extinction (for example, in the case where 0 < a1b2 < a2c1), the stochastic system (1.4) may be permanent. In fact, we shall see in the subsequent numerical simulations that that strong noise on the intra- specific interaction rate of the predator population may lead to the permanence of the prey, which is different from the deterministic case. Theorems 2.2 and 2.4 means that the signs of λ1 and λ2 determine the asymptotic behavior of the solution Xx(t) to (1.4) with initial value x ∈ R2,◦ + . A natural question is: what happens to the solution Xx(t) when λ1λ2 = 0? In fact, the set Ξ = {(a1, a2, b1, b2, c1, c2, µ1, µ2, ρ1, ρ2) ∈ R2×R8,◦ + : λ1λ2 = 0} has Lebesgue measure 0 in the space R2×R8,◦ + , and hence is negligible in the sense of the Lebesgue measure. Nevertheless, it is very interesting to explore the asymptotic behavior of the solution of (1.4) with parameters (a1, a2, b1, b2, c1, c2, µ1, µ2, ρ1, ρ2) ∈ Ξ. Unfortunately, this question remains open and we have to resort to new techniques. Finally, we do believe that our methods are applicable to stochastic predator-prey models with different types of functional response functions (i.e., replacing X1(t)X2(t) of system (1.4) by some function of X1(t) and X2(t)) (see [14, 15, 24, 31]) as well as to stochastic models with Markovian switching [12, 25]. 6 S. LI, S. GUO EJDE-2022/32 3. Numerical simulations To illustrate our theoretical results, we perform some numerical simulations of the solutions to (1.4), from which we find some interesting phenomena completely different from its corresponding deterministic system. By comparing the trajecto- ries of the stochastic system (1.4) with those of the associated deterministic system, we clarify the effect of the intensity of each while noise on the dynamics of system (1.4) and provide some reasonable biological interpretation. As we shall see in Theorems 2.2 and 2.4, and Figures 9 and 13, different from the associated deter- ministic predator-prey system, strong noise on the intra-specific interaction rate of prey results in the extinction of the predator species even though the predator population persists in the associated deterministic system, while strong noise on the intra-specific interaction rate of the predator population may lead to the per- manence of the prey even though the prey population dies out in the associated deterministic system. Furthermore, the population of the prey (or, predator) could decrease when there are strong noises on its corresponding intra-specific interaction rate, and strong noise on the inter-specific interaction rate of the prey population (or, the predator population) results in its extinction. By applying Milstein scheme in [7], we have the following discretization system of model (1.4), Xk+1 = Xk +Xk (a1 − b1Xk − c1Yk) ∆t+X2 k [ µ1ξ1,k √ ∆t+ 1 2 µ2 1(ξ2 1,k − 1)∆t ] +XkYk [ ρ1ξ3,k √ ∆t+ 1 2 ρ2 1(ξ2 3,k − 1)∆t ] , Yk+1 = Yk + Yk (a2 − b2Yk + c2Xk) ∆t+ Y 2 k [ µ2ξ2,k √ ∆t+ 1 2 µ2 2(ξ2 2,k − 1)∆t ] +XkYk [ ρ2ξ3,k √ ∆t+ 1 2 ρ2 2(ξ2 3,k − 1)∆t ] , where ∆t is the time increment and ξ1,k, ξ2,k and ξ3,k (k = 1, 2, 3, . . . ) are inde- pendent Gaussian random variables which follow the standard Normal distribution N(0, 1). We first illustrate Theorem 2.1 by the stochastic systems dXt = Xt(−0.1− 0.3Xt − 0.25Yt)dt+ 0.1X2 t dW1(t) + 0.1XtYtdW3(t), dYt = Yt(−0.1− 0.25Yt + 0.1Xt)dt+ 0.1Y 2 t dW2(t) + 0.1XtYtdW3(t), (3.1) and dXt = Xt(−0.1− 0.3Xt − 0.25Yt)dt+ µ1X 2 t dW1(t) + ρ1XtYtdW3(t), dYt = Yt(0.2− 0.25Yt + 0.1Xt)dt+ µ2Y 2 t dW2(t) + ρ2XtYtdW3(t). (3.2) It is easy to see that the deterministic system associated with (3.1) has a unique GAS equilibrium (0, 0), and that the deterministic system associated with (3.2) has a unique GAS equilibrium (0, 4/5). Obviously, Theorem 2.1(i) is illustrated by Figure 1, from which we see that both the prey and predator population go to extinction very fast and that the three white noises have no effect on sample paths of system (3.1). If the intensities of the three white noises is moderate, for example, µ1 = µ2 = ρ1 = ρ2 = 0.1, then the solution Xt of system (3.2) goes to zero while the solution Yt oscillates around the value 4 5 after some initial transients and the distribution of Yt converges weakly to π∗2 (see Figure 2). In fact, we can EJDE-2022/32 STOCHASTIC LOTKA-VOLTERRA PREDATOR-PREY MODELS 7 see that the larger the intensities of the white noises are, the larger the fluctuations of the solution Yt will be. However, too large intensity µ2 will make the solution Yt become more close to zero while too large intensity ρ1 will make the solution Xt vanish more rapidly (see Figures 3 and 4). (a) (b) (c) Figure 1. Numerical solutions to (3.1) with initial value (X0, Y0) = (0.4, 0.4): (a) trajectories of solutions; (b) phase por- trait; (c) top, right and inside of the box are the marginal and joint density distribution of solution (X,Y ), respectively. We next consider the stochastic system dXt = Xt(0.2− 0.3Xt − 0.25Yt)dt+ µ1X 2 t dW1(t) + ρ1XtYtdW3(t), dYt = Yt(−0.1− 0.25Yt + 0.1Xt)dt+ µ2Y 2 t dW2(t) + ρ2XtYtdW3(t), (3.3) whose associated deterministic system has a unique GAS equilibrium ( 2 3 , 0). We first consider system (3.3) with parameters µ1 = µ2 = ρ1 = ρ2 = 0.1. An easy cal- culation yields that λ1 ≈ −0.0363. Figure 5 shows that the solution of model (3.3) oscillates around the deterministic equilibrium ( 2 3 , 0) after some initial transients, that is, Yt goes to extinction and Xt weakly converges to π∗1 , which illustrates The- orem 2.2(i) and Corollary 2.3 as well. We also observe that increasing the intensity µ2 or ρ1 has little influence on the dynamical behaviours of (X,Y ) of (3.3). Fix the intensities µ2 = ρ1 = ρ2 = 0.1 and increase the intensity µ1 from 0.1 to 1, or fix the intensities µ1 = µ2 = ρ1 = 0.1 and increase the intensity ρ2 from 0.1 to 1, then we have λ1 ≈ −0.0647 and λ1 ≈ −0.2538, respectively. This, together with The- orem 2.2(i), implies that the predator population dies out. This theoretical result 8 S. LI, S. GUO EJDE-2022/32 (a) (b) (c) Figure 2. Numerical solutions to (3.2) with initial value (X0, Y0) = (0.4, 0.4): (a) trajectories of solutions; (b) phase por- trait; (c) top, right and inside of the box are the marginal and joint density distribution of solution (X,Y ), respectively. is illustrated by Figures 6 and 7, from which we can see that Yt goes to extinction faster if ρ2 becomes larger. Now we consider the stochastic system dXt = Xt(0.37− 0.3Xt − 0.25Yt)dt+ µ1X 2 t dW1(t) + ρ1XtYtdW3(t), dYt = Yt(−0.1− 0.25Yt + 0.24Xt)dt+ µ2Y 2 t dW2(t) + ρ2XtYtdW3(t), (3.4) whose associated deterministic system has a unique GAS equilibrium ( 47 54 , 98 225 ). We first consider system (3.4) with parameters µ1 = µ2 = ρ1 = ρ2 = 0.1. Direct calculation yields λ1 ≈ 0.1826. This, together with Theorem 2.2(ii) implies that the two species (X,Y ) coexist. This theoretical result is illustrated by Figure 8, from which we can see that the solution of model (3.4) oscillates around the deterministic equilibrium (47/54, 98/225) after some initial transients. We also observe that the increase of the intensity of µ2 and ρ1 has no much influence on the dynamical behaviours of solutions Xt and Yt of (3.4). Fix the intensities µ2 = ρ1 = ρ2 = 0.1 and increase µ1 from 0.1 to 2 (see Figure 9), or fix µ1 = µ2 = ρ1 = 0.1 and increase ρ2 from 0.1 to 1 (see Figure 10), then we have λ1 = −0.0277 and λ1 = −0.5552, respectively. This, together with Theorem 2.2(i), implies that Yt goes to extinction and Xt weakly converges to a boundary distribution π∗1 . Indeed, we can find out that the larger intensities µ1 and/or ρ2, the faster the convergence of Yt to 0, which means that strong noise whether on the intra-specific interaction rate of prey or EJDE-2022/32 STOCHASTIC LOTKA-VOLTERRA PREDATOR-PREY MODELS 9 (a) (b) (c) Figure 3. Numerical solutions to (3.2) with initial value (X0, Y0) = (0.4, 0.4) and parameters µ1 = ρ1 = ρ2 = 0.1 and µ2 = 1: (a) trajectories of solutions; (b) phase portrait; (c) top, right and inside of the box are the marginal and joint density dis- tribution of solution (X,Y ), respectively. inter-specific interaction rate of predator results in the extinction of the predator species. Now we consider the stochastic system dXt = Xt(0.2− 0.6Xt − 0.35Yt)dt+ µ1X 2 t dW1(t) + ρ1XtYtdW3(t), dYt = Yt(0.5− 0.5Yt + 0.3Xt)dt+ µ2Y 2 t dW2(t) + ρ2XtYtdW3(t), (3.5) whose associated deterministic system has a unique GAS equilibrium (0, 1). First, we choose the following parameters µ1 = µ2 = ρ1 = ρ2 = 0.1. An easy calculation yields λ1 ≈ 0.5 and λ2 ≈ −0.1515, which togethers with Theorem 2.4 implies that Xt goes to extinction and Yt converges to a boundary distribution π∗2 . From Figure 11 we can see that the solution of model (3.5) oscillates around the deterministic equilibrium (0, 1) after some initial transients. To find out the effect of intensity of noises on the dynamical behaviours of Xt and Yt, we shall increase the intensity µ1, µ2, ρ1 and ρ2, respectively. If we increase the intensities µ1 and ρ2, respectively, there is nothing new about the trajectories of Xt and Yt. If we fix µ1 = µ2 = ρ2 = 0.1 and increase intensity ρ1 from 0.1 to 1, we have λ1 ≈ 0.5 and λ2 ≈ −0.6416 and find that Xt converges to 0 faster (see Figure 12). This means that strong intensity of noise on the inter-specific interaction rate of the prey leads to a faster speed of extinction of the prey. Next, fix µ1 = ρ1 = ρ2 = 0.1 and increase µ2 from 0.1 to 10 S. LI, S. GUO EJDE-2022/32 (a) (b) (c) Figure 4. Numerical solutions to (3.2) with initial value (X0, Y0) = (0.4, 0.4) and parameters µ1 = µ2 = ρ2 = 0.1 and ρ1 = 2: (a) trajectories of solutions; (b) phase portrait; (c) top, right and inside of the box are the marginal and joint density dis- tribution of solution (X,Y ), respectively. 2, then we have λ2 ≈ 0.0818 and λ1 ≈ 0.5. Figure 13 shows that Xt is away from 0 and (Xt, Yt) coexist, which further illustrate Theorem 2.4(iii). This implies that strong noise on the intra-specific interaction rate of the predator population may lead to the permanence of the prey which is different from the deterministic case. Finally, we consider the system dXt = Xt(0.8− 0.6Xt − 0.35Yt)dt+ µ1X 2 t dW1(t) + ρ1XtYtdW3(t), dYt = Yt(0.5− 0.5Yt + 0.3Xt)dt+ µ2Y 2 t dW2(t) + ρ2XtYtdW3(t), (3.6) whose associated deterministic system has a unique GAS endemic equilibrium (5/9, 4/3). For the parameters µ1 = µ2 = ρ1 = ρ2 = 0.1, we have λ2 ≈ 0.4485 and λ1 ≈ 0.8868. This, together with Theorem 2.4(iii), implies that the two species (Xt, Yt) coexist. This theoretical result is illustrated by Figure 14, from which we can see that the solution of model (3.6) oscillates around the deterministic equilib- rium (5/9, 4/3). We next study the the effect of intensities of noises on the dynamical behaviours of (3.6). Fix µ2 = ρ1 = ρ2 = 0.1 and increase µ1 from 0.1 to 1.1, then we have λ2 ≈ 0.4485 and λ1 ≈ 0.7005. From Figure 15, we see that the two species coexist, but the solution Xt becomes closer to 0 and the probability distribution of the solution of Yt becomes closer to the invariant measure π∗2 . Fix µ1 = ρ1 = ρ2 = EJDE-2022/32 STOCHASTIC LOTKA-VOLTERRA PREDATOR-PREY MODELS 11 (a) (b) (c) Figure 5. Numerical solutions to (3.3) with initial value (X0, Y0) = (0.4, 0.4) and parameters µ1 = µ2 = ρ1 = ρ2 = 0.1: (a) Trajectories of solutions; (b) phase portrait; (c) top, right and inside of the box are the marginal and joint density distribution of solution (X,Y ), respectively. 0.1 and increase µ2 from 0.1 to 1, then we also observe that (X,Y ) coexists but the solution Xt is closer to the trajectories of ϕx(t) and the solution Yt is closer to 0 (see Figure 16). This implies that the population of the prey and predator could decrease when there are strong noises on their corresponding intra-specific interaction rate, respectively. Increasing the intensity ρ1 from 0.1 to 1 while fixing µ1 = µ2 = ρ2 = 0.1, we have λ1 ≈ 0.8868 and λ2 ≈ −0.0416, which satisfies the assumption of Theorem 2.4(i). In fact, we see from Figure 17 that the solution Xt is eventually extinct. As we increase the intensity ρ1 further, the solution Xt converges to 0 faster. In this case, strong noise on the inter-specific interaction rate of the prey population results in the extinction of the prey species. Similarly, increasing the intensity ρ2 from 0.1 to 1.1 while fixing µ1 = µ2 = ρ1 = 0.1, we have λ1 ≈ −0.1681 and λ2 ≈ 0.4485, which together with Theorem 2.4(ii) implies that the solution Yt goes to extinction. This theoretical result is illustrated by Figure 18. As we increase intensity ρ2 further, Yt converges to 0 faster. We see that strong noise on the inter-specific interaction rate of the predator population leads to the extinction of the predator population . 12 S. LI, S. GUO EJDE-2022/32 (a) (b) (c) Figure 6. Numerical solutions to (3.3) with initial value (X0, Y0) = (0.4, 0.4) and parameters µ1 = 1, µ2 = ρ1 = ρ2 = 0.1: (a) trajectories of solutions; (b) phase portrait; (c) top, right and inside of the box are the marginal and joint density distribution of solution (X,Y ), respectively. 4. Proof of Theorem 2.1 Lemma 4.1 ([12, 20, 21]). For every p ∈ (0, 3), there exists a positive constant Mp such that E ∫ t 0 ‖Xx(s)‖pds ≤Mp(t+ ‖x‖) for all x ∈ R2,◦ + and t ≥ 0. Note that E ∣∣∣ ∫ T 0 µjXj,x(t)dWj(t) +ρjX3−j,x(t)dW3(t) ∣∣∣2 = E ∫ T 0 (µ2 jX 2 j,x(t) +ρ2 jX 2 3−j,x(t))dt, then by Lemma 4.1 and Chebyshev’s inequality, we see that for any ς > 0, there exists a positive constant M̂ such that P {∣∣∣ ∫ T 0 µjXj,x(t)dWj(t) + ρjX3−j,x(t)dW3(t) ∣∣∣ ≤ M̂ ς √ T‖z‖ } ≥ 1− ς. (4.1) We define the stopping time of Xj,x as τσj,x = inf{t ≥ 0 : Xj,x(t) ≥ σ}, j = 1, 2. For R > 1 and δ > 0, let DR,δ 1 = (0, δ]× [R−1, R], DR,δ 2 = [R−1, R]× (0, δ]. Then we have the following results. Lemma 4.2. For any j = 1, 2, R > 1, T > 1, ς > 0 and σ > 0, there is a positive constant δ such that P{τσj,x ≥ T} ≥ 1− ς for all x ∈ DR,δ j . EJDE-2022/32 STOCHASTIC LOTKA-VOLTERRA PREDATOR-PREY MODELS 13 (a) (b) (c) Figure 7. Numerical solutions to (3.3) with initial value (X0, Y0) = (0.4, 0.4) and parameters µ1 = µ2 = ρ1 = 0.1, ρ2 = 1: (a) trajectories of solutions; (b) phase portrait; (c) top, right and inside of the box are the marginal and joint density distribution of solution (X,Y ), respectively. The proof of the above lemma is similar to that in [23] and hence we omit it. The following result means that Xj,x(t) is close to ψj,x(t) if X3−j,x(t) is small for a sufficiently long time. Lemma 4.3. For any j = 1, 2, R, T > 1 and positive constants ς, ν, and γ, there exists a positive constant σ such that for all x ∈ DR,σ 3−j, P{|ψj,x(t)−Xj,x(t)| < ν and |Φj,x(t)| < γ for all t ∈ [0, T ∧ τσ3−j,x] } ≥ 1− ς, where Φj,x(t) = 1 ψj,x(t) − 1 Xj,x(t) . The proof of the above lemma is similar to that in [23] and hence we omit it. Proof of Theorem 2.1(i). By exponential martingale inequality, we have P(Ωxj ) ≥ 1− ς, where Ωxj = {∫ t 0 ( µjXj,x(s)dWj(s) + ρjX3−j,x(s)dW3(s) ) ≤ ln 1 ς + 1 2 ∫ t 0 ( µ2 jX 2 j,x(s) + ρ2 jX 2 3−j,x(s) ) ds for all t ≥ 0 } , j = 1, 2. 14 S. LI, S. GUO EJDE-2022/32 (a) (b) (c) Figure 8. Numerical solutions to (3.4) with initial value (X0, Y0) = (0.4, 0.4) and parameters µ1 = µ2 = ρ1 = ρ2 = 0.1: (a) trajectories of X and Y ; (b) phase portrait; (c) top, right and inside of the box are the marginal and joint density distribution of solution (X,Y ), respectively. Applying Itô’s formula to the first equation of (1.4) yields that for ω ∈ Ωx1 and t > 0, lnX1,x(t) ≤ lnx+ ln 1 ς + a1t and so lim sup t→∞ 1 t lnX1,x(t) ≤ a1, (4.2) which implies that almost surely lim t→∞ X1,x(t) = 0, lim t→∞ 1 t ∫ t 0 X1,x(t)dt = 0. (4.3) Similarly, applying Itô’s formula to the second equation of (1.4) yields that for ω ∈ Ωx2 and t > 0, we have lnX2,x(t) ≤ lnx2 + ln 1 ς + a2t+ c2 ∫ t 0 X1,x(t)dt, (4.4) which together with (4.3) implies that almost surely limt→∞X2,x(t) = 0. The proof is complete. � Lemma 4.4. Assume that a1 < 0 < a2, then for any ς > 0, R > 1, there are T = T (ς, R) > 0 and δ0 = δ0(ς, R) > 0 such that P(Ω̌x) > 1 − ς for all x = (x1, x2) ∈ DR,δ0 2 , where Ω̌x = {lnX2,x(T ) − lnx2 ≥ 1 4a2T}. Moreover, there are EJDE-2022/32 STOCHASTIC LOTKA-VOLTERRA PREDATOR-PREY MODELS 15 (a) (b) (c) Figure 9. Numerical solutions to (3.4) with initial value (X0, Y0) = (0.4, 0.4) and parameters µ1 = 2, µ2 = ρ1 = ρ2 = 0.1: (a) trajectories of solutions; (b) phase portrait; (c) top, right and inside of the box are the marginal and joint density distribution of solution (X,Y ), respectively. T = T (ς) and δ1 = δ1(ς) > 0 such that lim sup N→∞ 1 N N−1∑ k=0 P{X2,x(kT ) ≤ δ1} ≤ ς for all x ∈ R2,◦ + . (4.5) Proof. As stated in the proof of Theorem 2.1(i), limt→∞X1,x(t) = 0 almost surely and hence for any given positive constant ς < a2/(2ρ 2 2) there exists T1 > 0 such that almost surely X1,x(t) < ς, 1 t ∫ t 0 X1,x(s)ds < ς, 1 t ∫ t 0 X2 1,x(s)ds < ς for all t > T1. Choose σ = σ(ς, R) > 0 such that b2σ + 1 2µ 2 2σ 2 ≤ a2/4. Lemma 4.2 implies that there exists δ0 = δ0(ς, R) > 0 such that P(Ωx3) ≥ 1−ς for all x ∈ DR,δ0 2 , where Ωx3 = {τσ2,x ≥ T1}. Let T2 > 16M̂2R/(ς2a2 2), then it follows from (4.1) that P(Ωx4) ≥ 1− ς, where Ωx4 = {∣∣∣ ∫ T2 0 µ2X2,x(t)dW2(t) + ρ2X1,x(t)dW3(t) ∣∣∣ ≤ 1 4 a2T2 } . Let T = T1 ∧ T2, then applying Itô’s formula to the second equation of (1.4) yields that lnX2,x(T )− lnx2 ≥ 1 4a2T for x ∈ DR,δ0 2 and ω ∈ Ω̌x = Ωx3 ∩Ωx4 . Consequently, 16 S. LI, S. GUO EJDE-2022/32 (a) (b) (c) Figure 10. Numerical solutions to (3.4) with initial value (X0, Y0) = (0.4, 0.4) and parameters ρ2 = 1, µ1 = µ2 = ρ1 = 0.1: (a) trajectories of solutions; (b) phase portrait; (c) top, right and inside of the box are the marginal and joint density distribution of solution (X,Y ), respectively. we obtain P(Ω̌x) > 1− ς. Similarly to [23, Proposition 3.2], we conclude that there are T = T (ς) and δ1 = δ1(ς) > 0 such that (4.5) holds, which implies that X2,x is away from zero for all x ∈ R2,◦ + . The proof is complete. � Proof of Theorem 2.1(ii)). As in the proof of Theorem 2.1(i), limt→∞X1,x(t) = 0 almost surely. In view of (4.2), there exists a positive random variable C1 such that X1,x(t) ≤ C1 exp{(a1 + ε)t} for t > 0 and sufficiently small ε. It is easy to see that there exists R > 1 such that lim sup t→∞ P{Xx(t) ∈ C} ≥ 1− ς, (4.6) where C , {R−1 ≤ x1 ∨ x2 ≤ R}. Using a similar argument as the proof of [23, Proposition 4.1], there exists δ > 0 such that P(Ωx5) ≥ 1− ς for all x ∈ DR,δ 1 , where Ωx5 = {|Φ2,x(t)| ≤ 1}. By Lemma 4.4, there exist T > 1 and δ1 > 0 such that lim sup N→∞ 1 N N−1∑ i=0 P{X2,x(iT ) < δ1} ≤ ς. (4.7) Note that C1 , C\(DR,δ 1 ∪ {(x1, x2) : x2 < δ1}) is compact and that X(t) is not recurrent in R2,◦ + (see [23, Proposition 4.1] for the related proof). Then EJDE-2022/32 STOCHASTIC LOTKA-VOLTERRA PREDATOR-PREY MODELS 17 (a) (b) (c) Figure 11. Numerical solutions to (3.5) with initial value (X0, Y0) = (0.4, 0.4) and parameters µ1 = µ2 = ρ1 = ρ2 = 0.1: (a) Trajectories of solutions; (b) phase portrait; (c) top, right and inside of the box are the marginal and joint density distribution of solution (X,Y ), respectively. the non-degeneracy of the diffusion implies that X(t) is transient and hence that limt→∞ P{Xx(t) ∈ C1} = 0. This together with (4.6) and (4.7) implies that lim sup N→∞ 1 N N−1∑ i=0 P { Xx(iT ) ∈ DR,δ 1 } ≥ 1− 2ς, and hence that there exists i0 such that P{Xx(i0T ) ∈ DR,δ 1 } ≥ 1 − 3ς, which together with Markov property implies that P {|Φ2,x(t)| < 1} ≥ (1− 3ς)2 ≥ 1− 6ς for all x ∈ R2,◦ + . In what follows, we shall show that X2,x(t) converges to ψ2,x(t). It follows from Itô’s formula that d[e2ρtΦ2 2,x(t)] = e2ρth(ψ2,x(t), X1,x(t), X2,x(t))dt + e2ρtg(ψ2,x(t), X1,x(t), X2,x(t))dW3(t) + 2ρe2ρtΦ2 2,x(t) where ρ = min{a2,−a1 − ς} and h(ψ, x, y) = 2 ( c2x− ρ2 2x 2 ) 1 y ( 1 ψ − 1 y ) + 2µ2 2(ψ − y) ( 1 ψ − 1 y ) + ρ2 2 x2 y2 − 2a2 ( 1 ψ − 1 y )2 , 18 S. LI, S. GUO EJDE-2022/32 (a) (b) (c) Figure 12. Numerical solutions to (3.5) with initial value (X0, Y0) = (0.4, 0.4) and parameters µ1 = µ2 = ρ2 = 0.1, ρ1 = 1: (a) trajectories of solutions; (b) phase portrait; (c) top, right and inside of the box are the marginal and joint density distribution of solution (X,Y ), respectively. g(ψ, x, y) = 2ρ2 x y ( 1 ψ − 1 y ) . For each γ > 0, we have P(Ωx6) ≥ 1− ς, where Ωx6 = {∫ t 0 g(ψ2,x(s), X1,x(s), X2,x(s))dW3(s) ≤ γ +mγ ∫ t 0 g2(ψ2,x(s), X1,x(s), X2,x(s))ds } and mγ = 2 γ ln 1 ς . For ω ∈ Ωx6 , e2ρtΦ2 2,x(t) ≤ γ + ∫ t 0 e2ρsh(ψ2,x(s), X1,x(s), X2,x(s))ds +mγ ∫ t 0 e4ρsg2(ψ2,x(s), X1,x(s), X2,x(s))ds. Thus for ω ∈ Ωx5 ∩ Ωx6 , Ψ̃x(t) = e2ρtΦ2 2,x(t) satisfies Ψ̃x(t) ≤ γ + 5C2 1ρ 2 2 ∫ t 0 exp{2(a1 + ε)s}Ψ̃x(s)ds+ C2 1 ( ρ2 2 + c22 a2 + 4mγρ 2 2 )∫ t 0 ds ψ2 2,x(s) EJDE-2022/32 STOCHASTIC LOTKA-VOLTERRA PREDATOR-PREY MODELS 19 (a) (b) (c) Figure 13. Numerical solutions to (3.5) with initial value (X0, Y0) = (0.4, 0.4) and parameters µ1 = ρ1 = ρ2 = 0.1, µ2 = 2: (a) trajectories of solutions; (b) phase portrait; (c) top, right and inside are the marginal and joint density distribution of solution (X,Y ), respectively. + 4C2 1mγρ 2 2 ∫ t 0 Ψ̃x(s) exp{(2(ρ+ a1 + ε)t}ds. In view of (2.2), there is T = T (ς, R) > 0 such that P(Ωx7) ≥ 1−ς for all x = (x1, x2) with x2 ∈ [R−1, R], where Ωx7 = {1 t ∫ t 0 ds ψ2 2,x(s) < 2Q2,−2, t > T } . Thus, for all ω ∈ Ωx5 ∩ Ωx6 ∩ Ωx7 and t > T , Ψ̃x(t) ≤ m1(t) + ∫ t 0 m2(s)Ψ̃x(s)ds, where m1(t) = γ + 2Q2,−2C 2 1 ( ρ2 2 + c22 a2 + 4mγρ 2 2 ) t, m2(t) = C2 1ρ 2 2(5 + 4mγ) exp{(2(ρ+ a1 + ε)s}. It follows from Gronwall’s inequality that Ψ̃x(t) ≤ m1(t) + ∫ t 0 m1(s)m2(s) exp {∫ t s m2(r)dr } ds 20 S. LI, S. GUO EJDE-2022/32 (a) (b) (c) Figure 14. Numerics of solutions to (3.6) with initial value (X0, Y0) = (0.4, 0.4) and parameters µ1 = µ2 = ρ1 = ρ2 = 0.1: (a) trajectories of X and Y ; (b) phase portrait; (c) top, right and inside of the box are the marginal and joint density distribution of solution (X,Y ), respectively. for all ω ∈ Ωx5 ∩ Ωx6 ∩ Ωx7 and t > T , and hence that almost surely lim t→∞ Φ2 2,x(t) ≤ lim t→∞ e−2ρt [ m1(t) + ∫ t 0 m1(s)m2(s) exp {∫ t s m2(r)dr } ds ] = 0. The proof is complete. � 5. Proof of Theorem 2.2 This section is devoted to the case where the prey population has a positive intrinsic growth rate while the predator population has a negative intrinsic growth rate, that is, a1 > 0 > a2. We have the following result. Lemma 5.1. Assume that a1 > 0 and λ1 < 0, then for any R > 1, ς, γ′ > 0, there exists δ̃ > 0 such that for all x ∈ DR,δ̃ 2 , P { lim sup t→∞ 1 t lnX2,x(t) < 0 and |Φ1,x(t)| < γ′ } ≥ 1− 3ς. Proof. Note that for each R > 1, there is R̄ = R(ε, R, T ) > 1 such that P{R̄−1 ≤ X1,x(t) ≤ R̄ for all t ∈ [0, T ]} ≥ 1 − ε if x ∈ [R−1, R] × [0, R]. By Lemma 4.3, for any R, T > 1, ς, γ′ > 0, there exists σ̃ > 0 such that P{|Φ1,x(t)| < γ′ for all t ∈ [t, T ∧ τ σ̃2,x]} ≥ 1 − ς for all x ∈ DR,σ̃ 2 . Set F (x, y) = −c2(x + y) + 1 2ρ 2 2(x − y)2. EJDE-2022/32 STOCHASTIC LOTKA-VOLTERRA PREDATOR-PREY MODELS 21 (a) (b) (c) Figure 15. Numerical solutions to (3.6) with initial value (X0, Y0) = (0.4, 0.4) and parameters µ1 = 1.1, µ2 = ρ1 = ρ2 = 0.1: (a) trajectories of solutions; (b) phase portrait; (c) top, right and inside of the box are the marginal and joint density distribution of solution (X,Y ), respectively. Choose λ ∈ (0,−λ1) and ϑ ∈ (λ2 ,− λ1 2 ). Since ∫∞ 0 F (φ, 0)f∗1 (φ)dφ = a2 − λ1 < ∞, we can find κ1, κ2 ∈ (0, 1) and d > 0 such that∫ ∞ κ1 F (φ, κ1)f∗1 (φ)dφ ≥ a2 − λ1 − ϑ = a2 + λ+ 3d,∫ ∞ κ−1 2 F (φ, κ1)f∗1 (φ)dφ ≤ d. There exists T2 = T2(ς, R) such that for all s ≥ 0, x ∈ [R−1, R]× [0, R] and t ≥ T2, P {1 t ∫ t 0 1{κ1≤ψ1,x(s)}F (ψ1,x(s), κ1)ds− (a2 + λ) ≥ 2d ≥ 1 t ∫ t 0 1{κ−1 2 ≤ψ1,x(s)}F (ψ1,x(s), κ1)ds } > 1− ς, and hence P(Ωx8) ≥ 1− ς, where Ωx8 = { a2 − 1 t ∫ t 0 1{κ1≤ψ1,x(s)≤κ−1 2 } F (ψ1,x(s), κ1)ds ≤ −λ for all t ≥ T2 } . 22 S. LI, S. GUO EJDE-2022/32 (a) (b) (c) Figure 16. Numerical solutions to (3.6) with initial value (X0, Y0) = (0.4, 0.4) and parameters µ2 = 1, µ1 = ρ1 = ρ2 = 0.1: (a) trajectories of solutions; (b) phase portrait; (c) top, right and inside of the box are the marginal and joint density distribution of solution (X,Y ), respectively. Hence, for ω ∈ Ωx8 ∩ {ϑx ≥ T2}, we have 1 t ∫ t 0 ( a2 + c2X1,x(s)− ρ2 2 2 X2 1,x(s) ) ds ≤ −λ (5.1) for all t ∈ [T2, ϑx], where ϑx = inf{t > 0 : |Φ1,x(t)| ≥ γ′ ∧ (κ1κ 2 2)}. In view of (5.1), for ω ∈ Ωx1 ∩ Ωx8 ∩ {ϑx ≥ T2}, we have lnX2,x(t) ≤ lnx2 + ln 1 ς − λt for all t ∈ [T2, ϑx], which together with Lemma 4.2 implies that there exists δ̃ = δ̃(ς, R) so small that ln δ̃ + ln 1 ς − λT2 < ln σ̃ and P(Ωx9) ≥ 1 − ς for all x ∈ DR,δ̃ 2 , where Ωx9 = {ζx , ϑx∧τ σ̃2,x ≥ T2}. Consequently, P(Ω̃x) ≥ 1−3ς, where Ω̃x = Ωx1∩Ωx8∩Ωx9 . As a result, for all x ∈ DR,δ̃ 2 and ω ∈ Ω̃x, we have lnX2,x(t ∧ τ σ̃2,x) ≤ lnx2 + ln 1 ς − λ(t ∧ τ σ̃2,x) < ln σ̃ for all t ≥ T2. and so t ∧ τ σ̃2,x < τ σ̃2,x for all t ≥ T2, x ∈ DR,δ̃ 2 and ω ∈ Ω̃x, which means that τ σ̃2,x = ϑx =∞ for all x ∈ DR,δ̃ 2 and ω ∈ Ω̃x. The proof is complete. � Lemma 5.2. If a1 > 0 > a2 and λ1 < 0, then for any R > 1, ς > 0, there exists δ > 0 and T = T (ς, R) > 0 such that P(Ω̆x) = 1 − ς for all x ∈ DR,δ 1 , where EJDE-2022/32 STOCHASTIC LOTKA-VOLTERRA PREDATOR-PREY MODELS 23 (a) (b) (c) Figure 17. Numerical solutions to (3.6) with initial value (X0, Y0) = (0.4, 0.4) and parameters ρ1 = 1, µ1 = µ2 = ρ2 = 0.1: (a) trajectories of solutions; (b) phase portrait; (c) top, right and inside of the box are the marginal and joint density distribution of solution (X,Y ), respectively. Ω̆x = {lnX1,x(T ) − lnx ≥ 1 4a1T}. Moreover, there are T = T (ς) and δ̃2(ς) > 0 such that lim sup N→∞ 1 N N−1∑ k=0 P { X1,x(kT ) ≤ δ̃2 } ≤ ς for all x ∈ R2,◦ + . (5.2) Proof. By the exponential martingale inequality, we have P(Ωx10) ≥ 1− ς, where Ωx10 = {∫ t 0 µ2ψ2,x(s)dW2(s) ≤ ln 1 ς + 1 2 ∫ t 0 µ2 2ψ 2 2,x(s)ds } . It follows from Itô’s formula that for ω ∈ Ωx10 and t > 0, lnψ2,x(t) ≤ lnx2+ln 1 ς +a2t, and hence lim sup t→∞ 1 t lnψ2,x(t) ≤ a2 < 0, which implies that limt→∞ ψ2,x(t) = 0 almost surely. There exists T3 > 0 such that P(Ωx11) ≥ 1− ς, where Ωx11 = { 1 T3 ∫ T3 0 ( a1 − c1(ψ2,x(t) + ν̂)− ρ2 1 2 (ψ2,x(t) + ν̂)2 ) dt ≥ 3 4 a1 } . It follows from Lemma 4.3 that we can choose σ̆ = σ̆(ς, R) > 0 such that b1σ̆ + 1 2µ 2 1σ̆ 2 ≤ a1/4 and P(Ωx12) ≥ 1− ς, where Ωx12 = {|ψ2,x(t)−X2,x(t)| < ν̂ for all t ∈ 24 S. LI, S. GUO EJDE-2022/32 (a) (b) (c) Figure 18. Numerical solutions to (3.6) with initial value (X0, Y0) = (0.4, 0.4) and parameters ρ2 = 1.1, µ1 = µ2 = ρ1 = 0.1: (a) trajectories of solutions; (b) phase portrait; (c) top, right and inside of the box are the marginal and joint density distribution of solution (X,Y ), respectively. [0, T3 ∧ τ σ̆1,x]}. Lemma 4.2 implies that there exists δ = δ(ς, R) such that P(Ωx13) ≥ 1− ς for all x ∈ DR,δ 1 , where Ωx13 = {τ σ̆1,x ≥ T3}. Let T4 > 16M̂2R/(ς2a2 1), then it follows that P(Ωx14) ≥ 1− ς, where Ωx14 = {∣∣ ∫ T5 0 µ1X1,x(t)dW1(t) + ρ1X2,x(t)dW3(t) ∣∣ ≤ a1 4 T4 } . Let T = T3 ∧ T4, then applying Itô’s formula to the first equation of (1.4) yields that lnX1,x(T )− lnx1 > 1 4a1T for x ∈ DR,δ 1 and ω ∈ Ω̆x = ∩14 i=10Ωxi . Consequently, P(Ω̆x) > 1 − ς. Similarly to [23, Proposition 3.2], we conclude that there are T = T (ς) > 0 and δ̃2 = δ̃2(ς) > 0 such that (5.2) holds, which implies that X1,x is away from zero for all x ∈ R2,◦ + . The proof is complete. � Proof of Theorem 2.2(i). It is easy to see that there exists R > 1 such that lim sup t→∞ P{Xx(t) ∈ C} ≥ 1− ς, (5.3) where C , {R−1 ≤ x1 ∨ x2 ≤ R}. By Lemma 5.1, there exists δ̃ > 0 such that P{lim sup t→∞ 1 t lnX2,x(t) ≤ −λ and |Φ1,x(t)| < γ′} ≥ 1− 3ς (5.4) EJDE-2022/32 STOCHASTIC LOTKA-VOLTERRA PREDATOR-PREY MODELS 25 for all x ∈ DR,δ̃ 2 . By Lemma 5.2, there exist T5 > 1 and δ̃2 > 0 such that lim sup N→∞ 1 N N−1∑ i=0 P{X1,x(iT5) < δ̃2} ≤ ς. (5.5) Using a similar arguments as in the proof of Theorem 2.1(ii), we can find i0 such that P{Xx(i0T5) ∈ DR,δ̃ 2 } ≥ 1−3ς, which together with (5.4) and Markov property implies that P { lim sup t→∞ 1 t lnX2,x(t) ≤ −λ } ≥ (1− 3ς)2 ≥ 1− 6ς. Namely, X2,x(t) goes to extinction. The convergence of X1,x(t) to π∗1 can be ob- tained by a similar method to Theorem 2.1(ii) as well. The proof is complete. � Lemma 5.3. If a1 > 0 > a2 and λ1 > 0 then then for any ς > 0, R > 1, there are T = T (ς, R) > 0 and δ0 = δ0(ς, R) such that P(Ω̂x) > 1 − 4ς for all x = (x1, x2) ∈ DR,δ0 2 , where Ω̂x = {lnX2,x(T ) − lnx2 ≥ 1 7λ1T}. Moreover, there are T = T (ς) and δ2 = δ2(ς) > 0 such that lim sup N→∞ 1 N N−1∑ k=0 P{X2,x(kT ) ≤ δ2} ≤ ς for all x ∈ R2,◦ + . (5.6) Proof. In view of the definition of λ1, we have a2 + ∫ ∞ 0 [ c2(u− ν)− ρ2 2 2 (u+ ν)2 ] f∗1 (u)du ≥ 2λ1 3 for sufficiently small ν. It follows from the ergodicity of ψ(t) that there is T = T (ς, R) > 1746M̂2R/(ς2λ2 1) such that P { a2 + 1 T ∫ T 0 [ c2 ( ψ1,R−1(s)− ν ) − ρ2 2 2 (ψ1,R(s) + ν) 2 ] ds ≥ λ1 2 } ≥ 1− ς. By the uniqueness of solution, ψ1,R−1(t) ≤ ψ1,x(t) ≤ ψ1,R(t) a.s. for all x ∈ [R−1, R]. Thus, we have P(Ωx15) ≥ 1− ς, where Ωx15 = { a2 + 1 T ∫ T 0 [ c2 (ψ1,x(s)− ν)− ρ2 2 2 (ψ1,x(s) + ν) 2 ] ds ≥ λ1 2 } . In view of Lemma 4.3, there is σ > 0 such that b2σ+ 1 2µ 2 2σ 2 ≤ λ1 3 and P(Ωx16) ≥ 1−ς, where Ωx16 = { |ψ1,x(t)−X1,x(t)| < ν for all t ∈ [0, T ∧ τσ2,x] } . It follows from Lemma 4.2 that there exists δ0 = δ0(ς, R) such that P(Ωx17) ≥ 1− ς for all x ∈ DR,δ0 2 , where Ωx17 = {τσ2,x ≥ T}. Using a similar argument as the proof of Lemma 4.4, we have P(Ωx18) ≥ 1− ς, where Ωx18 = {∣∣ ∫ T 0 µ2X2,x(s)dW2(s) + ρ2X1,x(s)dW3(s) ∣∣ ≤ λ1 42 T } . Applying Itô’s formula to the second equation of (1.4) yields that for x ∈ DR,δ0 2 and ω ∈ Ω̂x = ∩18 i=15Ωxi , lnX2,x(T )− lnx2 ≥ ∫ T 0 [ a2 + c2 ( ψ1,x(s)− ν ) − ρ2 2 2 (ψ1,x(s) + ν)2 ] dt− 5λ1 14 26 S. LI, S. GUO EJDE-2022/32 ≥ λ1 7 T. Consequently, we obtain P(Ω̂x) > 1 − 4ς. Similarly to [23, Proposition 3.2], we conclude that there are T = T (ς) and δ2 = δ2(ς) > 0 such that (5.6) holds, which implies that X2,x is away from zero for all x ∈ R2,◦ + . � Proof of Theorem 2.2(ii). Lemma 5.3 implies that X2,x is away from zero for all x ∈ R2,◦ + . If P {ω : limt→∞X1,x(t) = 0} > 0, then it follows from (4.4) that P ( lim sup t→∞ 1 t lnX2,x(t) < 0 ) > 0, which contradicts Lemma 5.3, and hence that almost surly lim supt→∞X1,x(t) > 0. In view of Lemma 4.1, there exists a constant C > 0 such that 1 t ∫ t 0 E(Xp 1,x(s) +Xp 2,x(s))ds ≤ C. According to in [4, Theorem 4.14], there exists a stationary distribution for Xx(t). This completes the proof. � 6. Proof of Theorem 2.4 This section is devoted to the case where both the predator population and the prey population have positive intrinsic growth rates, that is, a1 > 0 and a2 > 0. Lemma 6.1. Assume that a2 > 0 and λ2 < 0, then for any R > 1, ς, γ > 0, there exists δ̂ > 0 such that for all x ∈ DR,δ̂ 1 , P { lim sup t→∞ 1 t lnX1,x(t) ≤ 0 and |Φ2,x(t)| < γ } ≥ 1− 3ς. The proof of the above lemma is similar to that of Lemma 5.1 and hence we omit it. By Lemmas 4.3 and 5.3, using a similar method to Theorem 2.2, we can prove Theorem 2.4(i)(ii). Note that if λ1 > 0 and λ2 > 0, there exist T7 > 1 and δ̃2 > 0 such that lim sup N→∞ 1 N N−1∑ i=0 P{|X1,x(iT7)| ∧ |X2,x(iT7)| < δ̃2} ≤ 2ς. (6.1) The proof of Theorem 2.4(iii) is similar to the proof of Theorem 2.2(ii). Acknowledgements. This work was partially supported by the National Nat- ural Science Foundation of China (Grant Nos. 12071446 & 12101578), Natural Science Foundation of Hubei (Grant No. 2021CFB167), China Postdoctoral Sci- ence Foundation (Grant No. 2020M682507), the Post-doctoral Innovative Research Positions in Hubei Province, China (Grant No. 260789), and the Fundamental Re- search Funds for the Central Universities, China University of Geosciences (Wuhan) (Grant No. CUGST2). EJDE-2022/32 STOCHASTIC LOTKA-VOLTERRA PREDATOR-PREY MODELS 27 References [1] L. J. S. Allen; An Introduction to Stochastic Processes with Applications to Biology, Person Prentice Hall, 2003. [2] L. Arnold, W. Horsthemke, J. W. Stucki; The influence of external real and white noise on the Lotka-Volterra model, Biometrical Journal, 21(5) (1979), 451–471. [3] Y. T. Cai, C. C. Wang, D. Fan; Stability and bifurcation in a delayed predator-prey model with Holling-type IV response function and age structure, Electronic Journal of Differential Equations, 2021 (2021), no. 42, 1-16. [4] M. Chen; From Markov chains to non-equilibrium particle systems, World Scientific, 2004. [5] N. H. Dang, N. H. Du, T. V. Ton; Asymptotic behavior of predator-prey systems perturbed by white noise, Acta Applicandae Mathematicae, 115 (2011), 351–370. [6] J. Gao, S. Guo; Global dynamics and spatio-temporal patterns in a two-species chemotaxis system with two chemicals, Zeitschrift für angewandte Mathematik und Physik 72 (2021), 25. [7] D. J. Higham; An algorithmic introduction to numerical simulation of stochastic differential equations, SIAM Review, 43 (2001), 525–546. [8] C. Ji, D. Jiang; Dynamics of a stochastic density dependent predator-prey system with Beddington-DeAngelis functional response, Journal of Mathematical Analysis and Applica- tions, 381 (2011), 441–453. [9] D. Jiang, C. Ji, X. Li, D. O’Regan; Analysis of autonomous Lotka-Volterra competition sys- tems with random perturbation, Journal of Mathematical Analysis and Applications, 390(2) (2012), 582–595. [10] S. Li, S. Guo; Permanence of a stochastic prey-predator model with a general functional response, Mathematics and Computers in Simulation, 187 (2021), 308–336. [11] S. B. Li, Y. Xiao, Y. Dong; Diffusive predator-prey models with fear effect in spatially het- erogeneous environment, Electronic Journal of Differential Equations, 2021 (2021), no. 70, 1-31. [12] X. Li, D. Jiang, X. Mao; Population dynamical behavior of Lotka-Volterra system under regime switching, Journal of Computational and Applied Mathematics, 232(2) (2009), 427– 448. [13] C. Liu, S. Guo; Steady states of Lotka-Volterra competition models with nonlinear cross- diffusion, Journal of Differential Equations, 292(2021), 247-286. [14] Q. Liu, D. Q. Jiang, T. Hayat, B. Ahmad; Stationary distribution and extinction of a sto- chastic predator-prey model with additional food and nonlinear perturbation, Applied Math- ematics and Computation, 320 (2018), 226–239. [15] C. Lu, X. H. Ding; Periodic solutions and stationary distribution for a stochastic predator- prey system with impulsive perturbations, Applied Mathematics and Computation, 350 (2019), 313–322. [16] Q. Luo, X. Mao; Stochastic population dynamics under regime switching, Journal of Mathe- matical Analysis and Applications, 334(1) (2007), 69–84. [17] J. Lv, K. Wang; Asymptotic properties of a stochastic predator-prey system with Holling II functional response, Communications in Nonlinear Science and Numerical Simulation, 16 (2011), 4037–4048. [18] L. Ma, S. Guo; Positive solutions in the competitive Lotka-Volterra reaction-diffusion model with advection terms, Proceedings of the American Mathematical Society, 149(7) (2021), 3013-3019 [19] L. Ma, S. Guo; Bifurcation and stability of a two-species reaction-diffusion-advection com- petition model, Nonlinear Analysis: Real World Applications 59 (2021), 103241 [20] X. Mao, G. Marion, E. Renshaw; Environmental Brownian noise suppresses explosions in population dynamics, Stochastic Processes and their Applications, 97(1) (2002), 95–110. [21] X. Mao, S. Sabais, E. Renshaw; Asymptotic behavior of stochastic Lotka-Volterra model, Journal of Mathematical Analysis and Applications, 287 (2003), 141–156. [22] J. Murray; Mathematical Biology I: An Introduction, (3rd Edition), Springer-Verlag, 2002. [23] D. H. Nguyen, G. Yin; Coexistence and exclusion of stochastic competitive Lotka-Volterra models, Journal of Differential Equations, 262 (2017), 1192–1225. [24] K. Nosrati, M. Shafiee; Dynamic analysis of fractional-order singular Holling type-II predator-prey system, Applied Mathematics and Computation, 313 (2017), 159–179. 28 S. LI, S. GUO EJDE-2022/32 [25] M. Q. Ouyang, X. Y. Li; Permanence and asymptotical behavior of stochastic prey-predator system with Markovian switching, Applied Mathematics and Computation, 266 (2015), 539– 559. [26] H. Qiu, S. Guo, S. Li; Stability and bifurcation in a predator-prey system with prey-taxis, International Journal of Bifurcation and Chaos, 30 (2) (2020), 2050022 [27] R. Rudnicki; Long-time behaviour of a stochastic prey-predator model, Stochastic Processes and their Applications, 108 (2003), 93–107. [28] R. Rudnicki, K. Pichór; Influence of stochastic perturbation on prey-predator systems, Math- ematical Biosciences, 206 (2007), 108–119. [29] R. Rudnicki, K. Pichór, M. Tyran-Kaminska; Markov semigroups and their applications. In: Garbaczewski, P., Olkiewicz, R. (eds.) Dynamics of Dissipation. Lecture Notes in Physics, Springer, Berlin, 597 (2002), 215–238. [30] F. Vadillo; Comparing stochastic Lotka-Volterra predator-prey models, Applied Mathematics and Computation, 360 (2019), 181–189. [31] A. Yagi, T. V. Ton; Dynamic of a stochastic predator-prey population, Applied Mathematics and Computation, 218 (2011), 3100–3109. [32] Y. Zhu, M. Liu; Permanence and extinction in a stochastic service-resource mutualism model, Applied Mathematics Letters, 69 (2017), 1–7. Shangzhi Li Center for Mathematical Sciences, China University of Geosciences, Wuhan, Hubei 430074, China Email address: lishangzhi@cug.edu.cn Shangjiang Guo (corresponding author) School of Mathematics and Physics, China University of Geosciences, Wuhan, Hubei 430074, China Email address: guosj@cug.edu.cn 1. Introduction 2. Main results 3. Numerical simulations 4. Proof of Theorem ?? 5. Proof of Theorem ?? 6. Proof of Theorem ?? Acknowledgements References