Electronic Journal of Differential Equations, Vol. 2025 (2025), No. 32, pp. 1–19. ISSN: 1072-6691. URL: https://ejde.math.txstate.edu, https://ejde.math.unt.edu DOI: 10.58997/ejde.2025.32 DYNAMIC BEHAVIOR OF A STOCHASTIC PREDATOR-PREY MODEL WITH STAGE-STRUCTURE AND NONLINEAR PERTURBATION XIANG ZHANG, MING KANG, FENGJIE GENG Abstract. In this article, we propose and analyze stage-structured stochastic predator-prey model, where a nonlinear perturbation is considered. Firstly, we prove that the stochastic system has a unique global positive solution. And then we discuss the ergodic stationary distribution of the random system. In addition, we obtain sufficient conditions for the extinction of populations. Finally, numerical simulations verify our theoretical results and show that nonlinear perturbation has more practical significance than linear perturbation. 1. Introduction In recent years, the study of population ecology has attracted extensive attention. Popula- tion ecology mainly studies the dynamical behavior of populations and the relationship between populations and the environment, which is of great significance to the survival of various species in nature and the sustainable utilization of environmental resources. There are three kinds of relationships among species: predator-prey relationship, competition relationship, and reciprocity relationship. The predator-prey relationship mainly describes the interaction between predators and prey which plays an important role in the development of population dynamics. Over the past few decades, predator-prey systems have attracted a lot of scholars’ attention, and several achievements have been produced [2, 4, 5, 7, 9, 10, 12, 19, 24, 26, 28]. As we know, the functional response function describes the biological transfer differences be- tween different species and powerfully affects the dynamical properties of the models. There- fore, many scholars have proposed various predator-prey systems with different functional re- sponses, such as Lotka-Volterra [28], Holling I-IV [10, 26], Beddington-DeAngelis [2, 5, 7, 19], Crowley–Martin [4, 9], ratio-dependent [12, 24] and so on. In addition to the functional response function, the stage structure is another crucial element to investigate predator-prey interactions. In the real world, since the reproduction and survival rate of biological populations are usually dependent on age or stage, their lives can be divided into two stages: juvenile and adult. In recent years, some scholars have devoted themselves to studying the predator-prey models with stage structure for prey or predator [1, 6, 8, 22, 26]. Zhao et al. [26] studied a stochastic predator–prey system with stage structure for prey and obtained the sufficient criteria for the existence of stationary distribution and ergodicity. Bai and Xu [1] investigated a stochastic predator–prey system with stage structure for predator and constructed sufficient conditions for global asymptotic stability. Meanwhile, in the actual ecosystem, environmental noises are everywhere and have a certain impact on the population. Therefore, it is more reasonable to investigate the law of population for the predator-prey model by virtue of stochastic differential equation. Up to now, there have been different kinds of approaches to introduce random perturbations [14, 16, 23, 25, 26]. Most scholars introduce the linear stochastic perturbation into the deterministic system to reveal the influence 2020 Mathematics Subject Classification. 34F05, 34E10, 34D20, 92B05. Key words and phrases. Predator-prey model; stage structure; nonlinear perturbation; global positive solution; stationary distribution. ©2025. This work is licensed under a CC BY 4.0 license. Submitted September 7, 2024. Published April 2, 2025. 1 2 X. ZHANG, M. KANG, F. GENG EJDE-2025/32 of environmental noises [13, 23, 25, 26, 27]. Especially, Zhang et al. [25] recently proposed the following predator-prey model with different response functions to juvenile and adult prey: dx(t) = (ry −mx− αxz − d1x)dt+ σ1xdB1(t), dy(t) = (mx− sy2 − βyz (1 + ay)(1 + bz) − d2y)dt+ σ2ydB2(t), dz(t) = (pαxz + qβyz (1 + ay)(1 + bz) − δz2 − d3z)dt+ σ3zdB3(t), . (1.1) with initial value x(0) > 0, y(0) > 0, z(0) > 0, (1.2) where x(t), y(t), and z(t) denote the population densities of juvenile prey, adult prey and predator at time t, respectively. σ2 i (i = 1, 2, 3) denote the intensities of environmental noise and Bi(t) are the mutually independent standard Brownian motions. All the parameters are positive constants and their specific biological significance is shown in Table 1. In addition, p and q are constants, 0 < p, q < 1. For this model, they established sufficient conditions for the ergodic stationary distribution and extinction of the populations (1.1). Motivated by papers [11, 17, 18, 21], we adopt the way of nonlinear stochastic perturbation to describe the effects of more complicated noises on population dynamics. So far, there is less research on this aspect in the predator-prey models. Keeping this viewpoint in mind and reflecting the stage structure, we propose the following nonlinear stochastic differential equation according to system (1.1) and nonlinear perturbation theory: dx1(t) = (rx2 −mx1 − αx1y − d1x1)dt+ x1(σ11 + σ12x1)dB1(t), dx2(t) = (mx1 − sx2 2 − βx2y (1 + ax2)(1 + by) − d2x2)dt+ x2(σ21 + σ22x2)dB2(t), dy(t) = (pαx1y + qβx2y (1 + ax2)(1 + by) − δy2 − d3y)dt+ y(σ31 + σ32y)dB3(t), (1.3) with initial value x1(0) > 0, x2(0) > 0, y(0) > 0, where σ2 ij(i = 1, 2, 3; j = 1, 2, 3) denote the intensities of environmental noise, x1(t), x2(t), y(t) denote the population densities of juvenile prey, adult prey and predator at time t, respectively. The remaining parameters are the same as in system (1.1). Many scholars have paid attention to the impact of random perturbations on biological models, but most of them consider simple linear perturbations. With the increasingly complex living environment of organisms, such as human activities and global climate change, it is necessary and meaningful to study the effects of nonlinear perturbations on biological systems. The contents and methods of linear perturbation and nonlinear perturbation are similar, but it is more difficult to study nonlinear perturbation systems, for example, it is difficult to find suitable Lyapunov function and the inequalities used are more complex. It should be noted that the numerical simulation specifically displays that the secondary disturbance is more intense than the linear disturbance, which is closer to reality. The rest of this article is organized as follows. In Section 2, we demonstrate system (1.3) has a unique global positive solution which is a premise for the study of later questions. In Section 3, we obtain sufficient conditions for the existence and uniqueness of an ergodic stationary distribution. In Section 4, the sufficient conditions for the extinction of the prey and predator populations are established. In Section 5, the numerical simulations are provided to verify the derived theoretical results, meanwhile, the effect of high-order ambient noise is also revealed. 2. Existence and uniqueness of a global positive solution To study the properties of population dynamics, we first need to ensure the solution of system (1.3) is global and positive. EJDE-2025/32 DYNAMIC BEHAVIOR OF A STOCHASTIC PREDATOR-PREY MODEL 3 Theorem 2.1. For any given initial value (x1(0), x2(0), y(0)) ∈ R3 +, there is a unique solution (x1(t), x2(t), y(t)) ∈ R3 + of system (1.3) on t ≥ 0 and the solution will remain in R3 + with probability one, where R3 + = { (x1, x2, y) T ∈ R3 : x1 > 0, x2 > 0, y > 0 } . Proof. According to (1.3), since its coefficients satisfy the local Lipschitz condition in R3 +, then for any given initial value (x1(0), x2(0), y(0)) ∈ R3 + there is a unique local solution (x1(t), x2(t), y(t)) ∈ R3 + on t ∈ [0, τe), where τe indicates the explosion time. To prove the solution of system (1.3) is global, what we need to do is to show the explosion time τe = ∞ a.s.. Let n0 > 0 be sufficiently large such that initial value (x1(0), x2(0), y(0)) all lie within the interval [ 1 n0 , n0]. For each integer n ≥ n0, we define the stopping time τn = inf{t ∈ [0, τe) : min{x1(t), x2(t), y(t)} ≤ 1 n or max{x1(t), x2(t), y(t)} ≥ n}. We always set inf ∅ = ∞ (∅ denotes the empty set) in this paper. According to the definition of τn, it is easy to see τn is increasing as n → ∞. Denote τ∞ = limn→∞ τn a.s.. If τ∞ = ∞ a.s. is true, then τe = ∞ and (x1(t), x2(t), y(t)) ∈ R3 + a.s. for all t ≥ 0. Therefore, we just need to assure that τ∞ = ∞ a.s. is true, and then we can complete the proof. If the statement is false, then there exists a pair of constants T > 0 and ε ∈ (0, 1) such that P{τ∞ ≤ T} > ε. Thus, there exists an integer n1 ≥ n0, and for all n ≥ n1, such that P{τn ≤ T} ≥ ε. Define a C2-function V (x1, x2, y): R3 + → R+ as V (x1, x2, y) = 2 √ x1 − lnx1 + 2 √ x2 − lnx2 + 2 √ y − ln y. (2.1) Notice that 2 √ u− lnu ≥ 0, ∀u > 0. (2.2) Thus, for each (x1(t), x2(t), y(t)) ∈ R3 +, V (x1, x2, y) is a nonnegative function. Applying Itô’s formula [20] to V (x1, x2, y), one has dV (x1, x2, y) = LV (x1, x2, y)dt+ ( 1 √ x1 − 1 x1 ) (σ11 + σ12x1)x1dB1(t) + ( 1 √ x2 − 1 x2 ) (σ21 + σ22x2)x2dB2(t) + ( 1 √ y − 1 y ) (σ31 + σ32y)ydB3(t), where LV : R3 + → R is defined by LV = rx2 1 √ x1 −m √ x1 − α √ x1y − d1 √ x1 − √ x1(σ11 + σ12x1) 2 4 − r x2 x1 +m+ αy + d1 + (σ11 + σ12x1) 2 2 +mx1 1 √ x2 − sx 3 2 2 − β √ x2y (1 + ax2)(1 + by) − d2 √ x2 − √ x2(σ21 + σ22x2) 2 4 −m x1 x2 + sx2 + βy (1 + ax2)(1 + by) + d2 + (σ21 + σ22x2) 2 2 + pαx1 √ y + qβx2 √ y (1 + ax2)(1 + by) − δy 3 2 − d3 √ y − √ y(σ31 + σ32y) 2 4 − pαx1 − qβx2 (1 + ax2)(1 + by) + δy + d3 + (σ31 + σ32y) 2 2 ≤ rx2( 1 √ x1 − 1 x1 ) +mx1( 1 √ x2 − 1 x2 )− σ2 12x 5 2 1 4 + σ2 12x 2 1 + pα 2 x2 1 − σ2 22x 5 2 2 4 + σ2 22x 2 2 + sx2 + qβ 2 x2 2 − σ2 32y 5 2 4 + αy + pα 2 y + qβ 2 y + δy + σ2 32y 2 +m 4 X. ZHANG, M. KANG, F. GENG EJDE-2025/32 + d1 + d2 + d3 + β b + σ2 11 + σ2 21 + σ2 31 ≤ rx2 4 + mx1 4 − σ2 12x 5 2 1 4 + σ2 12x 2 1 + pα 2 x2 1 − σ2 22x 5 2 2 4 + σ2 22x 2 2 + sx2 + qβ 2 x2 2 − σ2 32y 5 2 4 + αy + pα 2 y + qβ 2 y + δy + σ2 32y 2 +m+ d1 + d2 + d3 + β b + σ2 11 + σ2 21 + σ2 31 = ( −σ2 12x 5 2 1 4 + σ2 12x 2 1 + pα 2 x2 1 + mx1 4 ) + ( −σ2 22x 5 2 2 4 + σ2 22x 2 2 + qβ 2 x2 2 + rx2 4 + sx2 ) + ( −σ2 32y 5 2 4 + σ2 32y 2 + αy + pα 2 y + qβ 2 y + δy ) +m+ d1 + d2 + d3 + β b + σ2 11 + σ2 21 + σ2 31 ≤ K1 +K2 +K3 +m+ d1 + d2 + d3 + β b + σ2 11 + σ2 21 + σ2 31, where K1 = sup x1∈R+ { −σ2 12x 5 2 1 4 + σ2 12x 2 1 + pα 2 x2 1 + mx1 4 } , K2 = sup x2∈R+ { −σ2 22x 5 2 2 4 + σ2 22x 2 2 + qβ 2 x2 2 + sx2 + rx2 4 } , K3 = sup y∈R+ { −σ2 32y 5 2 4 + σ2 32y 2 + αy + pα 2 y + qβ 2 y + δy } . Obviously, there exists a positive constant K satisfying LV (x, y, z) ≤ K. So dV (x1, x2, y) ≤ Kdt+ ( 1 √ x1 − 1 x1 ) (σ11 + σ12x1)x1dB1(t) + ( 1 √ x2 − 1 x2 ) (σ21 + σ22x2)x2dB2(t) + ( 1 √ y − 1 y ) (σ31 + σ32y)ydB3(t). (2.3) Integrating both sides of (2.3) from 0 to τn ∧ T and then taking the expectations on both sides leads to E [V (x1(τn ∧ T ), x2(τn ∧ T ), y(τn ∧ T ))] ≤ V (x1(0), x2(0), y(0)) +KT. (2.4) Note that for every ω ∈ {τn ≤ T} there is at least one of x1(τn, ω), x2(τn, ω), y(τn, ω) equal n or 1 n . Hence V (x1(τn, ω), x2(τn, ω), y(τn, ω)) ≥ min { 2 √ n− lnn, 2 √ 1 n + lnn } . Based on (2.4), we have V (x1(0), x2(0), z(0)) +KT ≥ E [ I{τn≤T}(ω)V (x1(τn, ω), x2(τn, ω), y(τn, ω)) ] ≥ εmin { 2 √ n− lnn, 2 √ 1 n + lnn } , where I{τn≤T}(ω) is the indicator function of {τn ≤ T}. Then letting n → +∞ results in the following contradiction +∞ > V (x1(0), x2(0), y(0)) +KT = +∞. Therefore we infer that τe = ∞ a.s. The proof is complete. □ EJDE-2025/32 DYNAMIC BEHAVIOR OF A STOCHASTIC PREDATOR-PREY MODEL 5 3. Existence of a stationary distribution It is known that random perturbations can destroy the stability of equilibrium states in deter- ministic systems, resulting in a stable distribution of random weak stability. In this section, we shall investigate the asymptotic behavior of the solution (x1(t), x2(t), y(t)) for system (1.3) and establish sufficient condition for its existence of a unique ergodic stationary distribution. Next, we introduce the following lemmas to state and prove the existence and uniqueness of the stationary distribution. Lemma 3.1 ([15]). The Markov process X(t) has a unique ergodic stationary distribution π(·) if there exists a bounded open domain U ⊂ Rd with regular boundary Γ, having the following properties: (1) the diffusion matrix A(x) is strictly positive definite for all x ∈ U ; (2) there exists a nonnegative C2-function V such that LV is negative for any Rd\U . Lemma 3.2 ([18]). For each x ≥ 0, we have x4 ≥ (3 4 x2 − 1 4 ) (x2 + 1). Theorem 3.3. Assume that 2s2 d2 ≥ σ2 22, δ2 d3 ≥ σ2 32 hold. Then (1.3) admits a unique stationary distribution π(·) and it is ergodic provided that λ > 0, where λ = 2 √ rm−m− d1 − β b − d2 − σ2 11 2 − σ2 21 − σ11σ12d2 + 2sm pαd2 ( d3 + σ2 31 + qβδ abd3 ) − 32r2σ2 12 27s2 − qαβ abd3 . Proof. We only need to show that conditions (1) and (2) in Lemma 3.1 hold. Firstly, we shall verify condition (1). The diffusion matrix of system (1.3) is A = (σ11 + σ12x1) 2 x2 1 0 0 0 (σ21 + σ22x2) 2 x2 2 0 0 0 (σ31 + σ32y) 2 y2  . It is apparent that A is positive definite for all (x1(t), x2(t), y(t)) ∈ R3 +, which means condition (1) in Lemma 3.1 holds. Next, we verify condition (2). Define a C2-function V : R3 + → R such that LV ≤ −1 on R3\U , where U is an open bounded set. Firstly, based on system (1.3), we know that L(− lnx1) = −r x2 x1 +m+ αy + d1 + (σ11 + σ12x1) 2 2 , (3.1) L(− lnx2) = −m x1 x2 + sx2 + βy (1 + ax2)(1 + by) + d2 + (σ21 + σ22x2) 2 2 ≤ −mx1 x2 + sx2 + β b + d2 + σ2 21 + σ2 22x 2 2. (3.2) We define V1(x1, x2) = − lnx1 − lnx2. 6 X. ZHANG, M. KANG, F. GENG EJDE-2025/32 Combining (3.1) and (3.2), we have LV1 = −r x2 x1 + αy +m+ d1 + (σ11 + σ12x1) 2 2 −m x1 x2 + sx2 + βy (1 + ax2)(1 + by) + d2 + (σ21 + σ22x2) 2 2 ≤ −rx2 x1 − mx1 x2 + αy +m+ d1 + σ2 11 2 + σ11σ12x1 + σ2 12 2 x2 1 + sx2 + β b + d2 + σ2 21 + σ2 22x 2 2 ≤ ( −2 √ rm+m+ d1 + β b + d2 + σ2 11 2 + σ2 21 ) + αy + σ11σ12x1 + σ2 12 2 x2 1 + sx2 + σ2 22x 2 2. (3.3) We define V2(x1) = u1(x1 + u2) v v , where u1 and u2 are positive constants which will be determined later, v ∈ (0, 1) is adequately small. Making use of Itô’s formula to function V2(x1) and according to Lemma 3.2, we obtain LV2 = u1(x1 + u2) v−1(rx2 −mx1 − αx1y − d1x1) − u1(1− v)(x1 + u2) v−2 x2 1(σ11 + σ12x1) 2 2 ≤ u1rx2 (x1 + u2) 1−v − u1(1− v)x2 1(σ11 + σ12x1) 2 2(x1 + u2) 2−v ≤ u1rx2 u1−v 2 − u1(1− v)σ2 12x 4 1 2(x1 + u2) 2−v ≤ u1rx2 u1−v 2 − u1u v+2 2 (1− v)σ2 12( x1 u2 ) 4 2(x1 u2 + 1) 2−v ≤ u1rx2 u1−v 2 − u1u v+2 2 (1− v)σ2 12( x1 u2 ) 4 4[(x1 u2 ) 2 + 1] ≤ u1rx2 u1−v 2 − u1u v+2 2 (1− v)σ2 12 4 [3 4 ( x1 u2 ) 2 − 1 4 ] = u1rx2 u1−v 2 − 3u1u v 2(1− v)σ2 12 16 x2 1 + u1u v+2 2 (1− v)σ2 12 16 , We choose u1 = 8 3(1−v)uv 2 , u2 = 8r 3(1−v)s , which yields LV2 ≤ sx2 + 32r2σ2 12 27(1− v) 2 s2 − σ2 12 2 x2 1. (3.4) We define V3(x1, x2, y) = V1 + V2 + 2s d2 y, EJDE-2025/32 DYNAMIC BEHAVIOR OF A STOCHASTIC PREDATOR-PREY MODEL 7 it then follows from (3.3), (3.4), and the condition of Theorem 3.3 that LV3 ≤ ( −2 √ rm+m+ d1 + β b + d2 + σ2 11 2 + σ2 21 + 32r2σ2 12 27(1− v) 2 s2 ) + αy + σ11σ12x1 + sx2 + σ2 22x 2 2 + sx2 + 2s d2 ( mx1 − sx2 2 − βx2y (1 + ax2)(1 + by) − d2x2 ) ≤ ( −2 √ rm+m+ d1 + β b + d2 + σ2 11 2 + σ2 21 + 32r2σ2 12 27(1− v) 2 s2 ) + αy + ( σ11σ12d2 + 2sm d2 )x1. (3.5) Let V4(x1, x2, y) = V3 + σ11σ12d2 + 2sm pαd2 (− ln y). Apply Itô’s formula to V4(x1, x2, y) and combine with (3.5), one obtains LV4 ≤ ( −2 √ rm+m+ d1 + β b + d2 + σ2 11 2 + σ2 21 + 32r2σ2 12 27(1− v) 2 s2 ) + αy + σ11σ12d2 + 2sm pαd2 ( − qβx2 (1 + ax2)(1 + by) + δy + d3 + (σ31 + σ32y) 2 2 ) + (σ11σ12d2 + 2sm d2 ) x1 + σ11σ12d2 + 2sm pαd2 (−pαx1) = ( −2 √ rm+m+ d1 + β b + d2 + σ2 11 2 + σ2 21 + 32r2σ2 12 27(1− v) 2 s2 ) + αy + σ11σ12d2 + 2sm pαd2 ( − qβx2 (1 + ax2)(1 + by) + δy + d3 + (σ31 + σ32y) 2 2 ) ≤ (−2 √ rm+m+ d1 + β b + d2 + σ2 11 2 + σ2 21 + σ11σ12d2 + 2sm pαd2 (d3 + σ2 31)) + 32r2σ2 12 27(1− v) 2 s2 + ( α+ σ11σ12d2 + 2sm pαd2 δ ) y + σ11σ12d2 + 2sm pαd2 σ2 32y 2. (3.6) We define V5(x1, x2, y) = V4 + ( α d3 + σ11σ12d2 + 2sm pαd2d3 δ)y. 8 X. ZHANG, M. KANG, F. GENG EJDE-2025/32 From (3.6) and the conditions of Theorem 3.3, it follows that LV5 ≤ −2 √ rm+m+ d1 + β b + d2 + σ2 11 2 + σ2 21 + σ11σ12d2 + 2sm pαd2 (d3 + σ2 31) + 32r2σ2 12 27(1− v) 2 s2 + ( α+ σ11σ12d2 + 2sm pαd2 δ ) y + σ11σ12d2 + 2sm pαd2 σ2 32y 2 + ( α d3 + σ11σ12d2 + 2sm pαd2d3 δ )( pαx1y + qβx2y (1 + ax2)(1 + by) − δy2 − d3y ) ≤ −2 √ rm+m+ d1 + β b + d2 + σ2 11 2 + σ2 21 + σ11σ12d2 + 2sm pαd2 (d3 + σ2 31) + 32r2σ2 12 27(1− v) 2 s2 + ( α d3 + σ11σ12d2 + 2sm pαd2d3 δ ) pαx1y + ( α d3 + σ11σ12d2 + 2sm pαd2d3 δ )qβ ab − ( α d3 + σ11σ12d2 + 2sm pαd2d3 δ ) δy2 + σ11σ12d2 + 2sm pαd2 σ2 32y 2 ≤ −2 √ rm+m+ d1 + β b + d2 + σ2 11 2 + σ2 21 + 32r2σ2 12 27(1− v) 2 s2 + qαβ abd3 + σ11σ12d2 + 2sm pαd2 ( d3 + σ2 31 + qβδ abd3 ) + ( α d3 + σ11σ12d2 + 2sm pαd2d3 δ ) pαx1y ≤ −λ(v) + (pα2 d3 + σ11σ12d2 + 2sm d2d3 δ ) x1y, (3.7) where λ(v) = 2 √ rm−m− d1 − β b − d2 − σ2 11 2 − σ2 21 − σ11σ12d2 + 2sm pαd2 ( d3 + σ2 31 + qβδ abd3 ) − 32r2σ2 12 27(1− v) 2 s2 − qαβ abd3 . Clearly, limv→0+ λ(v) = λ. By the continuity of the function λ(v) and λ > 0, we can pick v ∈ (0, 1) adequately small such that λ(v) > 0. Thus LV5 ≤ −λ+ (pα2 d3 + σ11σ12d2 + 2sm d2d3 δ ) x1y. We set V6(x1) = (σ11 + σ12x1) θ θ , V7(x2) = (σ21 + σ22x2) θ θ , V8(y) = 1 θ yθ, where 0 < θ < 1 is a sufficiently small constant. Then we have LV6 = rx2σ12(σ11 + σ12x1) θ−1 − (mx1 + αx1y + d1x1)σ12(σ11 + σ12x1) θ−1 − (1− θ)σ2 12 2 (σ11 + σ12x1) θ−2x2 1(σ11 + σ12x1) 2 ≤ rx2σ12(σ11 + σ12x1) θ−1 − (1− θ)σ2 12 2 (σ11 + σ12x1) θ−2x2 1(σ11 + σ12x1) 2 = rσ12x2 (σ11 + σ12x1) 1−θ − (1− θ)σ2 12x 2 1(σ11 + σ12x1) θ 2 ≤ rσ12x2 σ1−θ 11 − (1− θ)σθ+2 12 2 x2+θ 1 . (3.8) EJDE-2025/32 DYNAMIC BEHAVIOR OF A STOCHASTIC PREDATOR-PREY MODEL 9 In a similar way, LV7 ≤ mσ22x1 σ1−θ 21 − (1− θ)σθ+2 22 2 x2+θ 2 , (3.10) and LV8 = pαx1y θ + qβx2y θ (1 + ax2)(1 + by) − δy1+θ − d3y θ − (1− θ)(σ31 + σ32y) 2 2 yθ ≤ pαx1y θ + qβx2y θ − (1− θ)σ2 32 2 yθ+2. (3.9) We define a C2-function: V9(x1, x2, y) = MV5 + V6 + V7 + V8 − lnx2, where M > 0 satisfying −Mλ+B + d2 + σ2 21 ≤ −2, (3.12) where B = sup (x1,x2,y)∈R3 + { − (1− θ)σθ+2 12 4 x2+θ 1 − (1− θ)σθ+2 22 4 x2+θ 2 − (1− θ)σ2 32 4 y2+θ + σ11 θ−1σ12rx2 + σ21 θ−1σ22mx1 + pαx1y θ + qβx2y θ + sx2 + βy + σ2 22x 2 2 } . Furthermore, we note that V9(x1, x2, y) is not only continuous, but also tends to∞ as (x1(t), x2(t), y(t)) approaches the boundary of R3 +. So we can see that the function V9(x1, x2, y) has a minimum in the interior, which is denoted by V9(x 0 1, x 0 2, y 0). We define another C2-function V : R3 + → R by V (x1, x2, y) = V9(x1, x2, y)− V9(x 0 1, x 0 2, y 0). Based on inequalities (3.2), (3.6), (3.7), (3.8) and (3.10), one can obtain LV ≤ −Mλ+ (pα2 d3 + σ11σ12d2 + 2sm d2d3 δ ) Mx1y + σθ−1 11 σ12rx2 − (1− θ)σθ+2 12 2 x2+θ 1 + σθ−1 21 σ22mx1 − (1− θ)σθ+2 22 2 x2+θ 2 + pαx1y θ + qβx2y θ − (1− θ)σ2 32 2 y2+θ −m x1 x2 + sx2 + βy + d2 + σ2 21 + σ2 22x 2 2 ≤ −Mλ+ (pα2 d3 + σ11σ12d2 + 2sm d2d3 δ ) Mx1y − (1− θ)σθ+2 12 4 x2+θ 1 − (1− θ)σθ+2 22 4 x2+θ 2 − (1− θ)σ2 32 4 y2+θ −m x1 x2 +B + d2 + σ2 21 ≤ −2 + (pα2 d3 + σ11σ12d2 + 2sm d2d3 δ ) Mx1y − (1− θ)σθ+2 12 4 x2+θ 1 − (1− θ)σθ+2 22 4 x2+θ 2 − (1− θ)σ2 32 4 y2+θ −m x1 x2 . (3.10) Now we construct a bounded open domain Uε such that the condition (2) in Lemma 3.1 holds, namely, Uε = { (x1, x2, y) ∈ R3 + : ε < x1 < 1 ε , ε2 < x2 < 1 ε2 , ε < y < 1 ε } , where 0 < ε < 1 is a sufficiently small number. In the set R3 +\Uε, we can choose ε sufficiently small such that the following conditions hold:(pα2 d3 + σ11σ12d2 + 2sm d2d3 δ ) Mε ≤ 1, (3.11)(pα2 d3 + σ11σ12d2 + 2sm d2d3 δ ) Mε− (1− θ)σ2 32 4 ≤ 0, (3.12) 10 X. ZHANG, M. KANG, F. GENG EJDE-2025/32(pα2 d3 + σ11σ12d2 + 2sm d2d3 δ ) Mε− (1− θ)σ2 12 4 ≤ 0, (3.13) D − m ε ≤ −1, (3.14) D − (1− θ) 8 ( σ12 ε )θ+2 ≤ −1, (3.15) D − (1− θ) 4 ( σ22 ε2 )θ+2 ≤ −1, (3.16) D − (1− θ) 8 ( σ32 ε )θ+2 ≤ −1, (3.17) where D > 0 is defined later. For convenience, we can divide R3 +\Uε into six domains: U1 = { (x1, x2, y) ∈ R3 + : x1 ≤ ε } , U2 = { (x1, x2, y) ∈ R3 + : x1 > ε, x2 ≤ ε2 } , U3 = { (x1, x2, y) ∈ R3 + : y ≤ ε } , U4 = { (x1, x2, y) ∈ R3 + : x1 ≥ 1 ε } , U5 = { (x1, x2, y) ∈ R3 + : x2 ≥ 1 ε2 } , U6 = { (x1, x2, y) ∈ R3 + : y ≥ 1 ε } . Clearly, U c ε = R3 +\Uε = U1 ∪ U2 ∪ U3 ∪ U4 ∪ U5 ∪ U6. Next we shall validate LV ≤ −1 for each (x1, x2, y) ∈ U c ε . Case 1: For each (x1(t), x2(t), y(t)) ∈ U1, because x1y ≤ εy ≤ ε(1 + y2+θ) and inequalities (3.10), (3.11) and (3.12), we obtain LV ≤ −2 + (pα2 d3 + σ11σ12d2 + 2sm d2d3 δ ) Mx1y − (1− θ)σθ+2 12 4 x2+θ 1 − (1− θ)σθ+2 22 4 x2+θ 2 − (1− θ)σ2 32 4 y2+θ −m x1 x2 ≤ −2 + (pα2 d3 + σ11σ12d2 + 2sm d2d3 δ ) Mε(1 + y2+θ)− (1− θ)σ2 32 4 y2+θ ≤ −2 + 1 = −1. (3.18) Case 2: For each (x1(t), x2(t), y(t)) ∈ U2, considering (3.10) and (3.14), one can reach LV ≤ −2 + (pα2 d3 + σ11σ12d2 + 2sm d2d3 δ ) Mx1y − (1− θ)σθ+2 12 4 x2+θ 1 − (1− θ)σθ+2 22 4 x2+θ 2 − (1− θ)σ2 32 4 y2+θ −m x1 x2 ≤ −2 + (pα2 d3 + σ11σ12d2 + 2sm d2d3 δ ) Mx1y − (1− θ)σθ+2 12 8 x2+θ 1 − (1− θ)σ2 32 8 y2+θ − m ε ≤ D − m ε ≤ −1, (3.19) where D = sup (x1,x2,y)∈R3 + {(pα2 d3 + σ11σ12d2 + 2sm d2d3 δ ) Mx1y − (1− θ) 8 (σθ+2 12 x2+θ 1 + σ2 32y 2+θ) } . EJDE-2025/32 DYNAMIC BEHAVIOR OF A STOCHASTIC PREDATOR-PREY MODEL 11 Case 3: For each (x1(t), x2(t), y(t)) ∈ U3, similarly, by (3.10), (3.11) and (3.13), we have LV ≤ −Mλ+ (pα2 d3 + σ11σ12d2 + 2sm d2d3 δ ) Mx1y − (1− θ)σθ+2 12 4 x2+θ 1 − (1− θ)σθ+2 22 4 x2+θ 2 − (1− θ)σ2 32 4 y2+θ −m x1 x2 +B ≤ −2 + (pα2 d3 + σ11σ12d2 + 2sm d2d3 δ ) Mε(1 + x2+θ 1 )− (1− θ)σ2 12 4 x2+θ 1 = −2 + ((pα2 d3 + σ11σ12d2 + 2sm d2d3 δ ) Mε− (1− θ)σ2 12 4 ) x2+θ 1 + (pα2 d3 + σ11σ12d2 + 2sm d2d3 δ ) Mε ≤ −2 + 1 = −1. (3.20) Case 4: For each (x1(t), x2(t), y(t)) ∈ U4, by (3.10), (3.13) and (3.15), one has LV ≤ −2 + (pα2 d3 + σ11σ12d2 + 2sm d2d3 δBig)Mx1y − (1− θ)σθ+2 12 4 x2+θ 1 − (1− θ)σθ+2 22 4 x2+θ 2 − (1− θ)σ2 32 4 y2+θ −m x1 x2 ≤ D − (1− θ)σθ+2 12 8 (1 ε )2+θ ≤ −1. (3.21) Case 5: For any (x1(t), x2(t), y(t)) ∈ U5, by (3.10) and (3.16), one achieves LV ≤ −2 + (pα2 d3 + σ11σ12d2 + 2sm d2d3 δBig)Mx1y − (1− θ)σθ+2 12 4 x2+θ 1 − (1− θ)σθ+2 22 4 x2+θ 2 − (1− θ)σ2 32 4 y2+θ −m x1 x2 ≤ D − (1− θ)σθ+2 22 4 ( 1 ε2 )2+θ ≤ −1. (3.22) Case 6: For each (x1(t), x2(t), y(t)) ∈ U6, by (3.10) and (3.17), one obtains LV ≤ −2 + (pα2 d3 + σ11σ12d2 + 2sm d2d3 δ ) Mx1y − (1− θ)σθ+2 12 4 x2+θ 1 − (1− θ)σθ+2 22 4 x2+θ 2 − (1− θ)σ2 32 4 y2+θ −m x1 x2 ≤ D − (1− θ)σθ+2 32 8 (1 ε )2+θ ≤ −1. (3.23) Therefore, the condition A2 in Lemma 3.1 is satisfied. As a consequence, system (1.3) has a stationary distribution π(·) and it has the ergodic property. □ Remark 3.4. From the expression of λ, we may know that the nonlinear noise disturbance cannot be ignored. Only when σ12, σi1(i = 1, 2, 3) are small enough, the ergodic stationary distribution of system (1.3) may be established. The σ22 and σ32 seem to be irrelevant to the ergodic results. Their role in system (1.3) will be studied in the future. 12 X. ZHANG, M. KANG, F. GENG EJDE-2025/32 4. Extinction In this section, we will establish sufficient conditions for the extinction of all populations. Theorem 4.1. Let (x1(t), x2(t), y(t)) ∈ R3 + be the solution of (1.3) with any initial value (x1(0), x2(0), y(0)) ∈ R3 +. If min{d1 +m, d2}(R0 − 1)I{R0≤1} +max{d1 +m, d2}(R0 − 1)I{R0>1} − 1 2(σ−2 11 + σ−2 21 ) < 0, where R0 = √ rm (d1+m)d2 . Then the prey populations and predator population will die out, namely, lim t→∞ x1(t) = 0, lim t→∞ x2(t) = 0, lim t→∞ y(t) = 0. Proof. Let M0 = [ 0 r d1+m m d2 0 ] . Because M0 is a nonnegative irreducible matrix, it has a left eigenvector (ω1, ω2) corresponding to R0 [3], where (ω1, ω2) = (md2 , R0). Clearly, R0(ω1, ω2) = (ω1, ω2)M0. We define a C2-function V (x1, x2): R2 + → R by V (x1, x2) = k1x1 + k2x2, (4.1) where k1 = ω1 d1+m , k2 = ω2 d2 . Applying Ito’s formula to lnV (x1, x2), we obtain d(lnV ) = L(lnV )dt+ k1(σ11 + σ12x1)x1 V dB1(t) + k2(σ21 + σ22x2)x2 V dB2(t), where L(lnV ) = k1 V (rx2 −mx1 − αx1y − d1x1)− k21(σ11 + σ12x1) 2 x2 1 2V 2 + k2 V (mx1 − sx2 2 − βx2y (1 + ax2)(1 + by) − d2x2)− k22(σ21 + σ22x2) 2 x2 2 2V 2 = 1 V { ω1 d1 +m (rx2 −mx1 − αx1y − d1x1) } − k21(σ11 + σ12x1) 2 x2 1 2V 2 + 1 V {ω2 d2 (mx1 − sx2 2 − βx2y (1 + ax2)(1 + by) − d2x2) } − k22(σ21 + σ22x2) 2 x2 2 2V 2 . (4.2) By properties of inequalities, we have L(lnV ) ≤ 1 V { ω1r d1 +m x2 − ω1x1 + ω2m d2 x1 − ω2x2 } − k21σ 2 11x 2 1 + k21σ 2 12x 4 1 2V 2 − k22σ 2 21x 2 2 + k22σ 2 22x 4 2 2V 2 = 1 V (ω1, ω2)(M0(x1, x2) T − (x1, x2) T )− k21σ 2 11x 2 1 + k21σ 2 12x 4 1 2V 2 − k22σ 2 21x 2 2 + k22σ 2 22x 4 2 2V 2 = 1 V {(R0 − 1)ω1x1 + (R0 − 1)ω2x2} − k21σ 2 11x 2 1 + k21σ 2 12x 4 1 2V 2 − k22σ 2 21x 2 2 + k22σ 2 22x 4 2 2V 2 = (R0 − 1) V {k1(m+ d1)x1 + k2d2x2} − k21σ 2 11x 2 1 + k21σ 2 12x 4 1 2V 2 EJDE-2025/32 DYNAMIC BEHAVIOR OF A STOCHASTIC PREDATOR-PREY MODEL 13 − k22σ 2 21x 2 2 + k22σ 2 22x 4 2 2V 2 ≤ min {d1 +m, d2} (R0 − 1)I{R0≤1} +max {d1 +m, d2} (R0 − 1)I{R0>1} − k21σ 2 11x 2 1 + k21σ 2 12x 4 1 2V 2 − k22σ 2 21x 2 2 + k22σ 2 22x 4 2 2V 2 . Then, by the Cauchy inequality [20], we obtain V 2 = ( k1σ11x1 1 σ11 + k2σ21x2 1 σ21 )2 ≤ (k21σ 2 11x 2 1 + k22σ 2 21x 2 2) ( 1 σ2 11 + 1 σ2 21 ) . (4.3) Hence, d(lnV ) ≤ { min {d1 +m, d2} (R0 − 1)I{R0≤1} +max{d1 +m, d2}(R0 − 1)I{R0>1} − 1 2(σ−2 11 + σ−2 21 ) − k21σ 2 12x 4 1 2V 2 − k22σ 2 22x 4 2 2V 2 } dt+ k1(σ11 + σ12x1)x1 V dB1(t) + k2(σ21 + σ22x2)x2 V dB2(t). Integrating from 0 to t and then dividing by t on both sides, we obtain lnV (t)− lnV (0) t ≤ min {d1 +m, d2} (R0 − 1)I{R0≤1} − 1 2(σ−2 11 + σ−2 21 ) + max {d1 +m, d2} (R0 − 1)I{R0>1} − 1 t ∫ t 0 k21σ 2 12x 4 1 2V 2 ds − 1 t ∫ t 0 k22σ 2 22x 4 2 2V 2 ds+ 1 t ∫ t 0 k1σ11x1 V dB1(s) + 1 t ∫ t 0 k2σ21x2 V dB2(s) + 1 t ∫ t 0 k1σ12x 2 1 V dB1(s) + 1 t ∫ t 0 k2σ22x 2 2 V dB2(s). (4.4) For convenience, let M1(t) := ∫ t 0 k1σ11x1 V dB1(s), M2(t) := ∫ t 0 k2σ21x1 V dB2(s), M3(t) := ∫ t 0 k1σ12x 2 1 V dB1(s), M4(t) := ∫ t 0 k2σ22x 2 2 V dB2(s). By the strong law of large numbers for martingales [20], we have lim t→∞ M1(t) t = 0, lim t→∞ M2(t) t = 0 a.s. (4.5) In addition, applying exponential martingales inequality [20], it is easy to see that for any positive constants T , u and v, which has P { sup 0≤t≤T [ Mi(t)− u 2 ⟨Mi(t),Mi(t)⟩ ] > v } ≤ e−uv, i = 3, 4. (4.6) Choosing T = k, u = 1, v = 2 ln k, one obtains P { sup 0≤t≤k [ Mi(t)− 1 2 ⟨Mi(t),Mi(t)⟩ ] > 2 ln k } ≤ 1 k2 , i = 3, 4. (4.7) On the basis of the Borel-Cantelli lemma [20], for almost all ω ∈ Ω, there is a random integer k0 = k0(ω) such that for k ≥ k0, which yields sup 0≤t≤k [ Mi(t)− 1 2 ⟨Mi(t),Mi(t)⟩ ] ≤ 2 ln k, i = 3, 4. (4.8) 14 X. ZHANG, M. KANG, F. GENG EJDE-2025/32 That is, M3 ≤ 2 ln k + 1 2 ⟨M3(t),M3(t)⟩ = 2 ln k + 1 2 ∫ t 0 ( k1σ12x 2 1 V ) 2 ds, (4.9) M4 ≤ 2 ln k + 1 2 ⟨M4(t),M4(t)⟩ = 2 ln k + 1 2 ∫ t 0 ( k2σ22x 2 2 V ) 2 ds. (4.10) For all 0 ≤ t ≤ k, k ≥ k0, substituting (4.9) and (4.10) into (4.4), we obtain lnV (t)− lnV (0) t ≤ min {d1 +m, d2} (R0 − 1)I{R0≤1} − 1 2(σ−2 11 + σ−2 21 ) + max {d1 +m, d2} (R0 − 1)I{R0>1} + M1(t) t + M2(t) t + 4 ln k t . For 0 ≤ k − 1 ≤ t ≤ k, it follows that lnV (t)− lnV (0) t ≤ min{d1 +m, d2}(R0 − 1)I{R0≤1} − 1 2(σ−2 11 + σ−2 21 ) + 4 ln k k − 1 + max{d1 +m, d2}(R0 − 1)I{R0>1} + M1(t) t + M2(t) t . Taking the limit superior and based on (4.5), we obtain lim sup t→∞ lnV (t) t ≤ min {d1 +m, d2} (R0 − 1)I{R0≤1} − 1 2(σ−2 11 + σ−2 21 ) + max{d1 +m, d2}(R0 − 1)I{R0>1} < 0. (4.11) Therefore, lim sup t→∞ lnx1(t) t < 0, lim sup t→∞ lnx2(t) t < 0. It indicates that lim t→∞ x1(t) = 0, lim t→∞ x2(t) = 0. (4.12) Then there exists T0 > 0 and a set Ωε ⊂ R3 + such that P (Ωε) > 1 − ε and qβx2y (1+ax2)(1+by) ≤ εqβy a.s. for t ≥ t0 and ω ∈ Ωε. Using Ito’s formula, it follows that d(ln y(t)) = 1 y [ ραx1y + qβx2y (1 + ax2)(1 + by) − δy2 − d3y ] dt+ σ31dB3(t) + σ32ydB3(t)− (σ31 + σ32y) 2 2 dt ≤ ( ραε+ qβε− d3 − σ2 31 2 ) dt+ σ31dB3(t) + σ32ydB3(t)− σ2 32y 2 2 dt. (4.13) Integrating from 0 to t on both sides of (4.13) and then dividing by t we obtain ln y(t)− ln y(0) ≤ ( pαε+ qβε− d3 − σ2 31 2 ) t+ σ31B3(t)− ∫ t 0 σ2 32y 2 2 ds+ ∫ t 0 σ32ydB3(t) = (pαε+ qβε− d3 − σ2 31 2 )t+ σ31B3(t)− ∫ t 0 σ2 32y 2 2 ds+M5(t). (4.14) Similarly, on the basis of the exponential martingales inequality, for all 0 ≤ t ≤ n and n ≥ n0, we have M5(t) ≤ 2 lnn+ 1 2 ⟨M5(t),M5(t)⟩ = 2 lnn+ 1 2 σ2 32 ∫ t 0 y2ds. (4.15) Substituting (4.15) into (4.14), it turns that ln y(t)− ln y(0) ≤ (pαε+ qβε− d3 − σ2 31 2 )t+ σ31B3(t) + 2 lnn. (4.16) EJDE-2025/32 DYNAMIC BEHAVIOR OF A STOCHASTIC PREDATOR-PREY MODEL 15 For 0 ≤ n− 1 ≤ t ≤ n, dividing by t on both sides of (4.16), it follows that ln y(t)− ln y(0) t ≤ ραε+ qβε− d3 − σ2 31 2 + σ31B3(t) t + 2 lnn n− 1 . Taking the limit superior, lim sup t→∞ ln y(t) t ≤ pαε+ qβε− d3 − σ2 31 2 < 0, which yields limt→∞ y(t) = 0. This completes the proof. □ According to the sufficient conditions of population extinction in Theorem 4.1 and reference [6], nonlinear random disturbance has little effect on total population extinction, and linear white noise disturbance of juvenile and adult prey will accelerate the population extinction. 5. Numerical simulations In this section, we shall verify theoretical results and state the impact of nonlinear disturbance by using numerical simulations. For the numerical simulations, we adopt Milstein’s higher-order method to give numerical simulations. The discretization transformation of the stochastic system (1.3) is as follows: xj+1 1 = xj 1 + (rxj 2 −mxj 1 − αxj 1y − d1x j 1)∆t+ xj 1(σ11 + σ12x j 1) √ ∆tξ1,j + xj 1 2 ( σ2 11 + 3σ11σ12x j 1 + 2σ2 12(x j 1) 2 ) (ξ21,j − 1)∆t, xj+1 2 = xj 2 + ( mxj 1 − s(xj 2) 2 − βxj 2y j (1 + axj 2)(1 + byj) − d2x j 2 ) ∆t + xj 2(σ21 + σ22x j 2) √ ∆tξ2,j + xj 2 2 ( σ2 21 + 3σ21σ22x j 2 + 2σ2 22(x j 2) 2 ) (ξ22,j − 1)∆t, yj+1 = yj + ( pαxj 1z + qβxj 2y j (1 + axj 2)(1 + byj) − δ(yj) 2 − d3y j ) ∆t + yj(σ31 + σ32y j) √ ∆tξ3,j + yj 2 ( σ2 31 + 3σ31σ32y j + 2σ2 32(y j) 2 ) (ξ23,j − 1)∆t, where the time increment ∆t > 0, σ2 ij(i = 1, 2, 3; j = 1, 2, 3) are the intensities of the white noise, and the ξ2i,j denote mutually independent Gaussian random variables which follow the distribution N(0, 1). We choose initial values and other parameters as follows in Table 1. Table 1. Parameter values Param. Description Values r birth rate of juvenile prey 0.8 m mature conversion rate of juvenile prey 0.3 α attack coefficient of predators on juvenile prey 0.3 d1 death rate of juvenile prey 0.05 s intra-specific competition coefficient among predators 0.1 β capture rate of predator on adult prey 0.3 a handling time 4 b magnitude of mutual interference among predators 3 d2 death rate of adult prey 0.1 p conversion efficiency of predator capturing juvenile prey 0.4 q conversion efficiency of predator capturing adult prey 0.3 δ intra-specific interference coefficient among predators 0.15 d3 natural mortality rate of predator populations 0.05 16 X. ZHANG, M. KANG, F. GENG EJDE-2025/32 Example 5.1. To obtain the existence of stationary distribution for system (1.3) numerically, we choose time step ∆t = 0.001, σ2 i,j = 0.01(i, j = 1, 2, 3) and other parameter values see Table 1. By a direct computation, we gain λ = 2 √ rm−m− d1 − β b − d2 − σ2 11 2 − σ2 21 − σ11σ12d2 + 2sm pαd2 (d3 + σ2 31 + qβδ abd3 ) − 32r2σ2 12 27s2 − qαβ abd3 ≈ 0.014 > 0, and it is easy to see that 2s2 d2 ≥ σ2 22, δ2 d3 ≥ σ2 32. In Figure 1, the left column shows the paths of x1, x2, y of system (1.3) with an initial value (x1(0), x2(0), y(0)) = (2, 2, 1). The intensity of white noise as σ2 i,j = 0.01(i, j = 1, 2, 3) and other parameter values are given in Table 1. The red line in the left figure represents the solution for the corresponding deterministic system (1.3) and the blue line represents the solution for the nonlinear disturbed system (1.3). The right column shows the histogram of the population probability density function of x1, x2, y. That is to say, the conditions of Theorem 3.3 are satisfied. Therefore, (1.3) admits a unique ergodic stationary distribution π(·), as shown in Figure 1. Figure 1. Density functions and the paths of x1, x2, y for (1.3) with σ2 i,j = 0.01(i, j = 1, 2, 3). Example 5.2. We choose the initial value (x1(0), x2(0), y(0)) = (2, 2, 1) and time step ∆t = 0.001. To show the conclusion of Theorem 4.1, in Figure 2, we choose σ11 = σ21 = 2, σ12 = σ22 = σ31 = σ32 = 0.01 and other corresponding parameter values see Table 1. By calculations, we obtain (R0 − 1)[min{d1 +m, d2}(R0 − 1)I{R0≤1} +max{d1 +m, d2}(R0 − 1)I{R0>1}] − 1 2(σ−2 11 + σ−2 21 ) ≈ −0.083 < 0. Thus the conditions of Theorem 3.3 hold. Thus we can obtain that the prey populations and predator populations die out exponentially with probability one. EJDE-2025/32 DYNAMIC BEHAVIOR OF A STOCHASTIC PREDATOR-PREY MODEL 17 Figure 2. Paths of x1, x2, y for system (1.3) with σ11 = σ21 = 2, σ12 = σ22 = σ31 = σ32 = 0.01. Example 5.3. In this example, it is mainly shown that nonlinear perturbations are higher than linear perturbations and can even change the state of species’ existence and extinction. Similarly, we choose time step ∆t = 0.001. To clarify the difference between linear and nonlinear systems, let the initial value (x(0), y(0), z(0)) = (2, 2, 1), σ1 = σ2 = σ3 = 0.1 in system (1.1); the initial value (x1(0), x2(0), y(0)) = (2, 2, 1), σ11 = σ21 = σ31 = 0.1, σ12 = σ22 = σ32 = 1 in system (1.3) and other parameter values are the same, as shown in Table 1. In Figure 3, taking the juvenile prey population as an example, it can be directly seen that the juvenile prey population of the linear system (1.1) exists, while the juvenile population of the nonlinear system (1.3) is extinct. In addition, the solution of (1.3) oscillates more strongly than the solution of (1.1) and even changes the state of existence and extinction of the population. It can be concluded that nonlinear perturbation is more in line with today’s increasingly complex environment and has more practical significance. 6. Conclusion In this paper, we have studied the dynamical behavior of a stochastic predator-prey model with nonlinear perturbation. We first establish sufficient criteria for the existence and uniqueness of an ergodic stationary distribution of positive solutions to the random system by constructing the appropriate stochastic Lyapunov function. The existence of a stationary distribution means that all populations coexist and are randomly persistent over long periods. And then we obtain suffi- cient conditions for the extinction of the prey and predator populations. Finally, our theoretical results are verified by numerical simulations. It is worth mentioning that compared with the linear disturbance, the sufficient condition for the stationary distribution of the positive solution is more stringent when the quadratic random white noise disturbance term is added, which indicates that a small nonlinear disturbance will affect the existence state of the population. Nowadays, biolog- ical populations are facing more and more severe survival challenges, and a small environmental disturbance will have a great impact on the number and stability of the population. Therefore, nonlinear perturbation is more suitable for today’s environment. 18 X. ZHANG, M. KANG, F. GENG EJDE-2025/32 Figure 3. Paths of juvenile prey population for systems (1.1) and (1.3) with σi = 0.1(i = 1, 2, 3); σi1 = 0.1(i = 1, 2, 3), σi2 = 1(i = 1, 2, 3). References [1] Bai, H.; Xu R.; Global stability of a stochastic predator-prey model with stage-structure, Chinese Quarterly Journal of Mathematics. 32 (2017), no. 4, 425. [2] Beddington, J.; Mutual interference between parasites or predators and its effect on searching efficiency, The Journal of Animal Ecology. (1975), 331–340. [3] Berman, A.; Plemmons, R.; Nonnegative matrices in the mathematical sciences, SIAM, 1994. [4] Crowley, P.; Martin, E.; Functional responses and interference within and between year classes of a dragonfly population, Journal of the North American Benthological Society. (1989), no. 3, 211–221. [5] DeAngelis, D; Goldstein, R.; O’Neill, R.; A model for tropic interaction, Ecology. 56 (1975), no. 4, 881–892. [6] Dubey, B.; Kumar, A.; Dynamics of prey–predator model with stage structure in prey including maturation and gestation delays, Nonlinear Dynamics. 96 (2019), 2653–2679. [7] Feng, X.; Liu, X.; Sun, C.; and Jiang, Y.; Stability and hopf bifurcation of a modified leslie–gower predator–prey model with smith growth rate and b–d functional response, Chaos, Solitons & Fractals. 174 (2023), 113794. [8] Fu, S.; Zhang, L.; Hu, P.;Global behavior of solutions in a lotka–volterra predator–prey model with prey-stage structure, Nonlinear Analysis: Real World Applications. 14 (2013), no. 5, 2027–2045. [9] Guin, L.; Pal, P.; Alzahrani, J.; Ali, N.; Sarkar, K.; Djilali, S.; Zeb, A.; Khan, I.; Eldin, S.; Influence of allee effect on the spatiotemporal behavior of a diffusive predator–prey model with crowley–martin type response function, Scientific Reports. 13 (2023), no. 1, 4710. [10] Holling, C.; The functional response of predators to prey density and its role in mimicry and population regulation, The Memoirs of the Entomological Society of Canada. 97 (1965), no. S45, 5–60. [11] Jiang, D.; Zhou, B.; Han, B.; Ergodic stationary distribution and extinction of a n-species gilpin–ayala com- petition system with nonlinear random perturbations, Applied Mathematics Letters. 120 (2021), 107273. [12] Jost, C.; Arino, O.; Arditi, R.; About deterministic extinction in ratio-dependent predator–prey models, Bul- letin of Mathematical Biology. 61 (1999), no. 1, 19–32. [13] Kang, M.; Geng, F.; Zhao, M.; Dynamical behaviors of a stochastic predator-prey model with anti-predator behavior, Journal of Applied Analysis & Computation. 13 (2023), no. 3, 1209–1224. [14] Kang, M.; Zhang, X.; Geng, F.; Ma, Z.; Stationary distribution of a stochastic three species predator-prey model with anti-predator behavior, Journal of Applied Mathematics and Computing. 70 (2024), no. 2, 1365-1393. [15] Khasminskii, R.; Stochastic stability of differential equations, vol. 66, Springer Science & Business Media, 2011. [16] Li, Y.; Gao, H.; Existence, uniqueness and global asymptotic stability of positive solutions of a predator–prey system with holling ii functional response with random perturbation, Nonlinear Analysis: Theory, Methods & Applications. 68 (2008), no. 6, 1694–1705. [17] Liu, Q.; Jiang, D.; Stationary distribution and extinction of a stochastic sir model with nonlinear perturbation, Applied Mathematics Letters. 73 (2017), 8–15. EJDE-2025/32 DYNAMIC BEHAVIOR OF A STOCHASTIC PREDATOR-PREY MODEL 19 [18] Liu, Q.; Jiang, D.; Hayat, T.; Alsaedi, A.; Ahmad, B.; Dynamical behavior of a higher order stochastically perturbed siri epidemic model with relapse and media coverage, Chaos, Solitons & Fractals. 139 (2020), 110013. [19] Liu, S.; Beretta, E.; A stage-structured predator-prey model of beddington-deangelis type, SIAM Journal on Applied Mathematics. 66 (2006), no. 4, 1101–1129. [20] Mao, X.; Stochastic differential equations and applications, Elsevier, 2007. [21] Mu, X.; Jiang, D.; and Hayat, T.; Analysis on dynamical behavior of a stochastic phytoplankton-zooplankton model with nonlinear perturbation, Mathematical Methods in the Applied Sciences. 46 (2023), no. 5, 5505–5520. [22] Wei, F.; Fu, Q.; Hopf bifurcation and stability for predator–prey systems with beddington–deangelis type func- tional response and stage structure for prey incorporating refuge, Applied Mathematical Modelling. 40 (2016), no. 1, 126–134. [23] Xu, C.; Yu, Y.; and Ren, G.; Dynamic analysis of a stochastic predator–prey model with crowley–martin functional response, disease in predator, and saturation incidence, Journal of Computational and Nonlinear Dynamics. 15 (2020), no. 7, 071004. [24] Yu, T.; Wang, Q.; Zhai, S.; Exploration on dynamics in a ratio-dependent predator-prey bioeconomic model with time delay and additional food supply, Mathematical Biosciences and Engineering. 20 (2023), no. 8, 15094–15119. [25] Zhang, S.; Yuan, S.; and Zhang, T.; A predator-prey model with different response functions to juvenile and adult prey in deterministic and stochastic environments, Applied Mathematics and Computation. 413 (2022), 126598. [26] Zhao, X.; Zeng, Z.; Stationary distribution of a stochastic predator–prey system with stage structure for prey, Physica A: Statistical Mechanics and its Applications. 545 (2020), 123318. [27] Zhou, B.; Han, B.; Jiang, D.; Hayat, T.; Alsaedi, A.; Stationary distribution, extinction and probability density function of a stochastic vegetation–water model in arid ecosystems, Journal of Nonlinear Science. 32 (2022), no. 3, 30. [28] Zu, L.; Jiang, D.; O’Regan, D.; Hayat, T.; Ahmad, B.; Ergodic property of a lotka–volterra predator–prey model with white noise higher order perturbation under regime switching, Applied Mathematics and Computation. 330 (2018), 93–102. Xiang Zhang School of Science, China University of Geosciences, Xueyuan Road, 100083, Beijing, China Email address: xiangzhang hb@163.com Ming Kang School of Science, China University of Geosciences, Xueyuan Road, 100083, Beijing, China Email address: mingkang ah@163.com Fengjie Geng (corresponding autor) School of Science, China University of Geosciences, Xueyuan Road, 100083, Beijing, China Email address: gengfengjie hbu@163.com 1. Introduction 2. Existence and uniqueness of a global positive solution 3. Existence of a stationary distribution 4. Extinction 5. Numerical simulations 6. Conclusion References