EUROPEAN JOURNAL OF PURE AND APPLIED MATHEMATICS 2025, Vol. 18, Issue 2, Article Number 6097 ISSN 1307-5543 – ejpam.com Published by New York Business Global 1 Dynamics of Competitive Stochastic Predator-Prey2 Model with Holling Type II under Small Random3 Immigration4 Mogtaba Mohammed1, Jawdat Alebraheem 1,∗, Ismail M. Tayel1,5 Muhamad Hifzhudin Noor Aziz26 1 Department of Mathematics, College of Sciences,Majmaah University, Al-Majmaah 11952,7 Saudi Arabia8 2 Institute of Mathematical Sciences, Faculty of Science Universiti Malaya 50603,9 Kuala Lumpur, Malaysia10 11 Abstract. In this paper, we introduce a novel competitive stochastic predator-prey model by means of Holling type II with small random immigration. A Series of analytical results are con- ducted on the model such as the boundedness of the model’s solution, stochastic mean square stability, and the conditions of the persistence and extinction of the predator-prey population. All analytical results were validated with numerical simulations by choosing different sets of small random immigration. Conditions on the random immigration under which the system is stochasti- cally stable are specified. In addition, small immigration has been shown to have a great impact on persistence and extinction of the system, which can play an important role in species conservation, especially by controlling the conditions of stochastic small immigration. 2020 Mathematics Subject Classifications: 92-11, 60H10, 34D2012 Key Words and Phrases: Stochastic predator-prey model, Persistence, Extinction, Small ran-13 dom immigration14 15 1. Introduction16 One of the most interesting topics in biomathematics is the study of the dynamical17 behavior of the so-called predator-prey models. Most certainly, the initiative work in18 this direction is that of Lotka [1] and Voltera [2]. However, some of the flaws with this19 model are that the interaction between the prey and predator populations is instantaneous,20 which may not be the case in many natural ecosystems; moreover, the model presupposes21 that the prey population has infinitely many resources to consume. Lotka-Voltera model22 has been developed in its structure to display more realistic features of the dynamics23 ∗Corresponding author. DOI: https://doi.org/10.29020/nybg.ejpam.v18i2.6097 Email addresses: mogtaba.m@mu.edu.sa (M. Mohammed), j.alebraheem@mu.edu.sa (J. Alebraheem), i.tayel.edu.sa (I. M. Tayel) https://www.ejpam.com 1 Copyright: © 2025 The Author(s). (CC BY-NC 4.0) M. Mohamed et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6097 2 of 19 by many researchers, see, for example, [3–9]. Due to many questions that come from24 ecological points of view, up to date many researchers are still working to answer such25 questions, for some examples of predator-prey models include the Allee effect, the fear26 effect, immigration, prey refuge, and cannibalism see, for example, [10–13].27 One of the important and recent external elements that is considered to play a crucial28 role in the dynamics of predator-prey models is immigration. This factor could lead to29 stabilizing unstable models as in [14], and also helps in understanding the persistence and30 extinction of some models that were not clear to explain [15, 16]. In general, many of these31 elements are described by means of deterministic system of ordinary differential equations,32 which have it is own shortcomings in considering natural fluctuations that could affect33 the ecological environment. This constraint has increased interest in stochastic modeling34 tools [17–21], which capture the natural unpredictability that affects ecological systems.35 Small random immigration has mostly impacted simpler predator-prey interactions, such36 as Holling type I models [22]. These studies often focus on specific elements, such as37 stochastic square stability, without completely addressing wider dynamical behaviors like38 persistence and extinction. This indicates a significant research gap: understanding how39 random immigration effects more complicated predator-prey systems, particularly those40 regulated by Holling type II responses, which better depict predator saturation and are41 ecologically more realistic.42 Our goal in this paper is to fill this gap by investigating a stochastic predator-prey43 model with Holling type II dynamics and extra ecological interactions, emphasizing how44 random immigration influences the boundedness of solutions, stochastic stability, persis-45 tence, and extinction of both prey and predator populations. To provide a more complete46 understanding of these complex dynamics, we use both rigorous mathematical analysis47 and numerical simulation.48 2. Problem Statement49 We consider a deterministic competitive predator-prey model with Holling type II50 which has been used in earlier works for various perspectives [23–26]. The model equations51 are defined below52 du dt = ρu ( 1− u k ) − αuv 1 + βu (1) 53 dv dt = −σv + γuv 1 + βu − λv2 (2) 54 u0 = u(0), v0 = v(0) (3) where u and v are the prey and predator populations, respectively. Recently, Jawdat et55 al. [26], studied system (1)-(3) with small immigration in the prey population. The model’s56 stability and coexistence and extinction conditions between prey and predator populations57 were investigated. However, the effects of small immigration for both populations were58 not included. Hence, in this paper, we aim to study the stochastic counterpart of system59 M. Mohamed et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6097 3 of 19 (1)-(3), more precisely, we consider small random immigration on both prey and predator60 populations by means of the following system:61  du = [ ρu ( 1− u k ) − αuv 1 + βu ] dt+ δudB1(t) dv = [ −σv + γuv 1 + βu − λv2 ] dt+ ϵvdB2(t) u0 = u(0), v0 = v(0), (4) (5) (6) In general, the mathematical model (4)-(6) represents prey-predator interactions, where62 the term ρu ( 1− u k ) is the logistic growth rate for the prey population, αuv 1+βu is the func-63 tional response (which presents the prey consumption rate by the predator population),64 the term −σv is the loss rate for the predator population, the term γuv 1+βu is the response65 of predator population density through consumption of the prey population, λv2 is the66 effect of intraspesific competition for predator population. δudB1(t) and ϵvdB2(t) refer to67 low-intensity stochastic processes that represent small random immigration in the popu-68 lations (prey and predator). For more detailed interpretation on the parameters, we have69 the following:70 • ρ is the rate of growth in the prey population.71 • k is the carrying capacity of the model.72 • α is the rate of catch prey population by the predator.73 • β is the handling rate of the predator.74 • σ is the natural death rate in the predator population.75 • γ is the effectiveness of transforming ate up prey into predator birth.76 • λ is the intraspesific competition between predators.77 • δdB1 and ϵdB2 are the small random immigration in the prey and predator popula-78 tions, respectively.79 The stochastic processes B1(t) and B2(t) are Brownian motions with the following prop-80 erties:81 • B1(0) = B2(0) = 0.82 • B1(t) and B2(t) are continuous functions with probability 1, for all t ∈ [0, T ].83 • B1(t) and B2(t) have stationary independent increments.84 • B1(t+ s)−B1(t) and B2(t+ s)−B2(t) are normally distributed with zero mean and85 variance t.86 M. Mohamed et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6097 4 of 19 3. Model Boundedness87 As it is one of main validations of the model, in this section, we use some probabilistic88 and analytic in inequalities to obtain the boundedness of the model. System (4)-(6) could89 be written as:90 dX = M(t,X)dt+N(t,X)dB X0 = X(0) (7) (8) such that91 M : [ [0, T ]× R+ ]2 → R N : [ [0, T ]× R+ ]2 → R M(t,X) and N(t,X) are given by M(t,X) = [ ρu ( 1− u k ) − αuv 1+βu −σv + γuv 1+βu − λv2 ] , N(t,X) = [ δu ϵv ] X = [ u v ] dB(t) = [ dB1(t) dB2(t) ] Lemma 1. For some positive constant C, the functions M(t,X) and N(t,X) in (7) satisfy92 the following linear growth condition93 E∥M(t,X)∥2 + E∥N(t,X)∥2 ≤ C(1 + E∥X∥2) (9) Proof. From (4) and the fact that ρ, k, α and β are positive, one sees that94 u ≤ u0 + ρ ∫ t 0 uds+ δ ∫ t 0 udB1(s). Taking the expectation on both sides and using Itô’s isometry, we have95 E|u| ≤ u0 + ρE ∫ t 0 |u|ds. From Grownall’s inequality, one gets96 E|u| ≤ u0e ρt (10) Now,97 E∥M(t,X)∥2 = E∥ρu ( 1− u k ) − αuv 1 + βu ∥2 + E∥ − σv + γuv 1 + βu − λv2∥2 ≤ ρ2E∥u∥2 + γ2E∥u∥2 + ∥v∥2 ≤ CE∥u∥2(1 + ∥v∥2) (11) M. Mohamed et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6097 5 of 19 From (10) and (11) we have98 E∥M(t,X)∥2 ≤ C(1 + ∥v∥2) ≤ C(1 + ∥u∥2 + ∥v∥2) ≤ C(1 + ∥X∥2). (12) It is easy to see that99 E∥N(t,X)∥2 = δ2E∥u∥2 + ϵ2E∥v∥2 ≤ C(1 + ∥u∥2 + ∥v∥2) ≤ C(1 + ∥X∥2). (13) From (12) and (13) we conclude the proof.100 Theorem 1. For all p > 1 and some C > 0, the following condition is true101 E sup 0≤t≤T ∥X∥p ≤ C. (14) Proof. From (7), we have102 X = X0 + ∫ t 0 M(t,X)ds+ ∫ t 0 N(t,X)dB(t), (15) it is known that see for example [27] for αk ∈ R, r ∈ [0,∞) and 1 ≤ p < ∞ , that103 r∑ i=1 αp i ≤ ( r∑ i=1 αi )p ≤ rp−1 r∑ i=1 αp i . (16) From (16), the sup over [0, T ] and the expectation we have the following:104 E sup 0≤t≤T ∥X∥q ≤3q−1 ( |X0|q + E sup 0≤t≤T ∣∣∣∣ ∫ t 0 M(t,X)dt ∣∣∣∣q + E sup 0≤t≤T ∣∣∣∣ ∫ t 0 N(t,X)dB(t) ∣∣∣∣q). (17) Using Hölder inequality for 1 q + 1 q q−1 = 1, we get105 E sup 0≤t≤T ∣∣∣∣ ∫ t 0 M(t,X)dt ∣∣∣∣q ≤ E ∣∣∣∣ ∫ T 0 M(t,X)dt ∣∣∣∣q ≤ E (∣∣∣∣ ∫ T 0 M(t,X)dt ∣∣∣∣q)q× 1 q × (∫ T 0 (1) q q−1 dt )q× q−1 q ≤ T q−1E (∣∣∣∣ ∫ T 0 M(t,X)dt ∣∣∣∣q) . (18) M. Mohamed et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6097 6 of 19 From Lemma 1 and (16) we have106 E ∫ T 0 ( |M(t,X)|2 ) q 2 dt ≤ EC q 2 ∫ T 0 ( 1 + ∥X∥2 ) q 2 dt ≤ C q 2 2 q 2 −1E ∫ T 0 (1 + ∥X∥q)dt ≤ C q 2 2 q 2 −1 ( T + E ∫ T 0 ∥X∥qdt ) . (19) From the inequality (18) and (19), we have107 E sup 0≤t≤T ∣∣∣∣ ∫ t 0 M(t,X)dt ∣∣∣∣q ≤ [C q 2 2 q 2 −1T q + C q 2 2 q 2 −1T q−1 ] E ∫ T 0 ∥X∥qdt. (20) To deal with the stochastic term, we employ Burkholder-Davis-Gundyand Hölder inequal-108 ities and it yields109 E sup 0≤t≤T ∣∣∣∣ ∫ t 0 N(t,X)dB(s) ∣∣∣∣q ≤ CBE (∫ T 0 |N(t,X)|2dt ) q 2 CB ( E ∫ T 0 |N(t,X|qdt )(∫ T 0 (1) q q−2dt ) q 2 × q−2 q ≤ CBT q−2 2 E ∫ T 0 |N(t,X)|qdt. (21) Again we use Lemma 1 and inequality (16) for p = q 2 to get110 E sup 0≤t≤T ∣∣∣∣ ∫ t 0 N(s,X)ds ∣∣∣∣ ≤ CBT q−2 2 C q 2 2 q 2 −1E ∫ T 0 (1 + ∥X∥q)dt ≤ CBT q−2 2 C q 2 2 q 2 −1 [ T + E ∫ T 0 ∥X∥qdt ] ≤ CBT q 2C q 2 2 q 2 −1 + CBT q−2 2 C q 2 2 q 2 −1E ∫ T 0 ∥X∥qdt. (22) From (17), (21) and (22) one uses Grownall’s inequality to complete the proof.111 4. Stochastic mean square stability112 Here, we investigate the stochastic mean square stability for the system (4)-(6) around113 the interior equilibrium point of the deterministic version (1)-(3). Note that our system114 has three non-negative equilibrium points:115 1. P0 = (0, 0)116 M. Mohamed et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6097 7 of 19 2. P1 = (k, 0)117 3. P ∗ = (u∗, v∗), where (u∗, v∗) are the real positive roots of the following system:118  ρu∗ ( 1− u∗ k ) − αu∗v∗ 1 + βu∗ = 0, −σv∗ + γu∗v∗ 1 + βu∗ − λv∗ 2 = 0. (23) (24) The stochastic mean square stability will be studied around the interior positive equilib-119 rium point (u∗, v∗), we have120  du = ( ρu ( 1− u k ) − −αuv 1 + βu ) dt+ δ(u− u∗)dB1, dv = ( −σv + γuv 1 + βu ) dt− λv2 + ϵ(v − v∗)dW2. u0 = u(0), v0 = v(0). (25) (26) (27) By setting x = u− u∗ and y = v − v∗, we transform the model into:121 dX = f(t,X)dt+ g(t,X)dB, (28) where122 X = [ x y ] , 123 f(t,X) = [ u∗ ∂k1(u ∗,u∗) ∂u + k1 u∗ ∂k1(u ∗,u∗) ∂v v∗ ∂k2(u ∗,v∗) ∂u v∗∂k2(u∗,v∗) ∂v + k2 ] [ x y ] , where124 k1 = ρ ( 1− u k ) − αv 1 + βu , k2 = −σ + γu 1 + βu λv. This gives125 f(t,X) = u∗ ( −ρ k + αβv∗ (1+βu)2 ) u∗ ( −α 1+βu∗ ) v∗ ( (1+βu∗)γ−γβu∗ (1+βu∗)2 ) −λv∗ [x y ] , and126 g(t,X) = [ δx 0 0 ϵy ] [ x y ] . For strictly positive t0, we define the cylindrical domain D = R2 × [t0,∞). Hence, we127 introduce the function φ ∈ C2(D) as φ : D → R+. As in [28] the generator of Markov128 M. Mohamed et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6097 8 of 19 process for equation (28) is given by129 Lφ(t,X) = ∂φ(t,X) ∂t + fT (t,X) ∂φ(t,X) ∂X + 1 2 Trace [ gT (t,X) ∂2φ(t,X) ∂X2 g(t,X) ] . (29) The following theorem states the mean-square stability for our system (28)130 Theorem 2. Assume that131 ρu∗ k > αβu∗v∗ (1 + βu∗)2 + δ2 2 , λv∗ > ϵ2 2 , (30) (31) then, system (23)-(27) is mean-square asymptotically stable.132 Proof. Let133 φ(t,X) = 1 2 ( a1x 2 + a2y 2 ) , be the Markovian process generator, where ai(i = 1, 2) are some parameters to be deter-134 mined later. Then we have the following135 ∂φ ∂t = 0, ∂φ ∂X = (a1x, a2y), (32) 136 f t(t,X) ∂φ ∂X = a1u 2 ( −ρ k + αβv∗ (1 +Bu∗)2 ) − a1αu ∗ 1 + βu∗ uv + a2αv ∗ (1 + βu∗)2 uv − a2λv ∗v∗, (33) and137 1 2 Trace [ gT (t,X) ∂2φ(t,X) ∂X2 t g(t,X) ] = 1 2 a1u 2 + 1 2 a2v 2ϵ2. (34) Substituting (32)-(34) into (29), we get138 L(φ(t,X)) = u∗ ( −ρ k + αβv∗ (1 +Bu∗)2 ) u2a1 − αu∗ 1 + βu∗ uva1 + αv∗ (1 + βu∗)2 uva2 − λv∗v∗a2 + 1 2 a1δ 2u2 + 1 2 a2ϵ 2v2 = ( −ρ k u∗ − αβu∗v∗ (1 +Bu∗)2 − 1 2 δ2 ) a1u 2 − ( λv∗ − 1 2 ϵ2 ) a2v 2 − ( αu∗ 1 + βu∗ a1 − αv∗ (1 + βu∗)2 a2 ) uv. (35) M. Mohamed et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6097 9 of 19 Let αu∗a1 = γv∗ 1+βu∗a2, this implies139 L(φ(t,X)) = ( −ρu∗ k − αβu∗v∗ (1 +Bu∗)2 − 1 2 δ2 ) a1u 2 − ( λv∗ − 1 2 ϵ2 ) a2v 2. (36) For L(φ(t,X)) < 0, then140 −ρu∗ k > αβu∗v∗ (1 + βu∗)2 + 1 2 δ2, λv∗ > 1 2 ϵ2. This completes the proof.141 Corollary 1. When conditions (30) and (31) are not satisfied, then the system (4)-(6) is142 unstable.143 Remark 1. The main factors that affect the system’s stability are the rate of growth in the144 prey population, the carrying capacity of the model, the rate of catching prey population by145 the predator, the handling rate of the predator, intraspecific competition between predators,146 and the small random immigration in the prey and predator populations.147 5. Persistence and extinction148 In this section, we will study the persistence for both prey and predator populations149 and the extinction for the prey population.150 5.1. Persistence of the prey population151 The following lemma was collected from [29].152 Lemma 2. lim sup t→∞ 0≤t≤T 1 t lnu ≤ 0 a.s.. (37) Theorem 3. The prey population u admits weak persistence in the average with probability153 almost if154 ρ− δ2 2 > 0. Proof. To prove this, we need to show that ul > a > 0, where155 ul := lim sup t→∞ 1 t ∫ t 0 uds, (38) we prove this by contradiction. Assume that ϱ1 be small enough such that156 −σ − 1 2 ϵ2 + γϱ1 < 0, (39) M. Mohamed et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6097 10 of 19 and157 ρ− 1 2 δ2 − ρ k ϱ1 > 0, for each ϱ1 > 0, there exists a solution (û, v̂) where158 P{ul < ϱ1}. (40) Therefore159 d ln v̂ ≤ [−σ − ϵ2 2 + γû]dt+ ϵdB2. Once again, we apply the integration over (0, t) and divide by t, we get160 ln v̂ − ln v̂0 t ≤ 1 t ∫ t 0 ( −σ − ϵ2 2 ) ds+ 1 t ∫ t 0 γûds+ 1 t ∫ t 0 ϵdB2(s), = −σ − ϵ2 2 + γ t 1 t ∫ t 0 γûds+ 1 t ∫ t 0 ϵdB2(s). (41) As before, from the strong law of large number, we get161 lim sup t→∞ ln v̂ t ≤ −σ − ϵ2 2 + γϱ1 < 0, (42) this implies that lim v̂ = 0. Therefore,162 d ln(û) = [ ρ− ρ k û− αûv̂ 1 + βû − 1 2 δ2 ] dt+ δdB1(t). (43) Integration followed by division by t, give163 ln û− ln û0 t = 1 t ∫ t 0 ( ρ− 1 2 δ2 ) ds− ρ kt ∫ t 0 ûds− α t ∫ t 0 ûv̂ 1 + βû ds+ δ t ∫ t 0 dB1(s), = ρ− 1 2 δ2 − ρ kt ∫ t 0 ûds− α t ∫ t 0 ûv̂ 1 + βû ds+ δ t ∫ t 0 dB1(s). Taking the lim supt→∞ in both sides with (39), (40) and the law of large number, we164 obtain165 lim sup 1 t ln û = ρ− δ2 2 − ρ k ϱ1 > 0. (44) This contradicts the lemma (2). Therefore, ul > 0. This completes the proof.166 5.2. Extinction of the prey population167 Theorem 4. The prey population u goes to extinct with probability almost surely when168 ρ− δ2 2 ≤ 0. (45) M. Mohamed et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6097 11 of 19 Proof. Let us construct the following comparison system;169 dU ≤ U ( ρ− ρU k ) dt+ δUdB1(t), U0 = u0. From Itô’s Lemma for the function φ(t, U) = lnU , we get170 d lnU = [ ρ− ρU k − 1 2 δ2 ] dt+ δdB1(t). (46) Applying integration on the interval (0, t) followed by division by t to get171 lnU − lnu0 t = 1 t ∫ t 0 ( ρ− ρu k − 1 2 δ2 ) ds+ 1 t ∫ t 0 δdB1(s). (47) Let us note that from the strong law of large number, one gets172 lim sup t→∞ 1 t ∫ t 0 δdB1(s) = 0, (48) and this implies that173 lim sup t→∞ lnU t ≤ ρ− 1 2 δ2 < 0 a.s. (49) Using the comparison theorem for SDEs, we have174 lim sup t→∞ 1 t lnu < 0, thus175 lim t→∞ u = 0. 5.3. Extinction for the predator population176 Theorem 5. If177 γk ( ρ− δ2 2 ) < ρ ( δ + ϵ2 2 ) , (50) then the predator population goes to extinct with probability almost surely.178 Proof. Let us consider the following two cases179 (i) ρ− δ2 2 < 0. It is clear that this will imply ul < 0. As in (41), we have180 ln v − ln v0 t ≤ −σ − ϵ2 2 + γ t ∫ t 0 uds+ ϵB2(t) t . (51) From this and the fact that ul < 0, then181 lim sup t→∞ ln v t ≤ −σ − ϵ2 2 < 0, then limt→∞v = 0, which means the extinction of predator population V (t).182 M. Mohamed et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6097 12 of 19 (ii) ρ − δ2 2 > 0. In this case, one can always ϱ2 > 0 at a time T1 such that for t > T1183 and has B1(t) t < ϱ2. Thus, we have184 lnu− lnu0 ≤ ∫ t 0 ( ρ− δ2 2 ) ds− ρ k ∫ t 0 uds+ δB1(t). ≤ ( ρ− δ2 2 + ϱ2 ) t− ρ k ∫ t 0 uds. From [30, Lemma 2.2], one obtains ul ≤ k ( ρ− δ2 2 + ϱ2 ) ρ Now, since ϱ2 → 0 as t → ∞, then we have ul ≤ k ( ρ− δ2 2 ) ρ . From (17), one sees that185 lim sup t→∞ 1 t ln v ≤ −σ − ϵ2 2 + γul ≤ −σ − ϵ2 2 + k ( ρ− δ2 2 ) ρ < 0. Thus limt→∞ v = 0.186 Corollary 2. The persistence of the predator population is obtained by reversing the con-187 dition (50).188 6. Numerical Simulations and discussion189 To validate our theoretical results presented in the previous sections, we performed several numerical simulations in this section. We employ the “StochasticRungeKuttaS- calarNoise” command in MATHEMATICA 11.1 to perform numerical simulations, as de- scribed on the Wolfram website. The initial populations values are chosen to be u0 = 0.9, v0 = 0.8. And the set of parameters are chosen such that the stability and coexistence between prey190 and predator populations hold for the deterministic model (1)-(3) in order to check the191 effect of the random immigration of the system as follows:192 ρ = 0.4, k = 2.0, α = 0.75, β = 0.5, σ = 0.2, γ = 0.6, λ = 0.1. (52) M. Mohamed et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6097 13 of 19 Since the model equations in (4)-(6) are taken to be dimensionless, the parameters were193 freely chosen. With these selections, our persistence equilibrium point is (u∗, v∗) ≈194 (0.5254, 0.4965). Now, to study the impact of small random immigration on the prey-195 predator populations (u, v), several considerations of δ and ϵ were taken into account196 (δ = 0, ϵ = 0), (δ = 0.05, ϵ = 0.05), (δ = 0.4, ϵ = 0.4), (δ = 0.4, ϵ = 1.2) and197 (δ = 1.0, ϵ = 0.1).198 199 0 20 40 60 80 100 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 Time P o p u la ti o n Prey Predator (a) Time series 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 Prey P re d a to r (b) Phase plane Figure 1: Dynamical behaviour for the deterministic model (1)-(3) with initial values (u0 = 0.9, v0 = 0.8.) and parameters values are given in (52). Figure1 represents the dynamical behavior of the deterministic model (1)-(3) at the200 equilibrium point (u∗, v∗) ≈ (0.5254, 0.4965), in which the system is stable and persist.201 For more detailed analysis of this model, we refer to [25].202 203 0 20 40 60 80 100 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 Time P o p u la ti o n Prey Predator (a) Time series 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 Prey P re d a to r (b) Phase plane Figure 2: Dynamical behaviour for the stochastic model (4)-(6) with initial values (u0 = 0.9, v0 = 0.8.), (δ = 0.05, ϵ = 0.05) and parameters values are given in (52). Figure 2 represents the stochastic model (4)-(6) with (δ = 0.05, ϵ = 0.05), in which the204 calculations give the following:205 • ρu∗ k − αβu∗v∗ (1+Bu∗)2 − δ2 2 ≈ 0.0425 > 0.206 M. Mohamed et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6097 14 of 19 • λv∗ − ϵ2 2 ≈ 0.0484 > 0.207 • ρ− δ2 2 ≈ 0.3988 > 0.208 • γk ( ρ− δ2 2 ) − ρ ( δ + ϵ2 2 ) ≈ 0.4478 > 0.209 Thus, the conditions of Theorems 2 and 3 are satisfied, which means the system (4)-(6)210 is stochastically stable and the prey population u is weakly persistence in average with211 probability almost surely and the predator population v achieved the persistence through212 simulation under the same condition.213 214 0 20 40 60 80 100 0 1 2 3 4 Time P o p u la ti o n Prey Predator (a) Time series 0.0 0.5 1.0 1.5 2.0 2.5 3.0 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Prey P re d a to r (b) Phase plane Figure 3: Dynamical behaviour for the stochastic model (4)-(6) with initial values (u0 = 0.9, v0 = 0.8.), (δ = 0.4, ϵ = 0.4) and parameters values are given in (52). Figure 3 indicates the stochastic model (4)-(6) with (δ = 0.4, ϵ = 0.4), in which the215 calculations give the following:216 • ρu∗ k − αβu∗v∗ (1+Bu∗)2 − δ2 2 ≈ −0.0363 < 0.217 • λv∗ − ϵ2 2 ≈ 0− 0.0305 < 0.218 • ρ− δ2 2 ≈ 0.32 > 0.219 • γk ( ρ− δ2 2 ) − ρ ( δ + ϵ2 2 ) ≈ 0.096 > 0.220 From the above calculations, one can see that the conditions of Theorem 2 are not met,221 and therefore the system is unstable. The condition of Theorem 3 is satisfied, therefore222 the prey population is weakly persistent in average with probability almost surely. In223 the meanwhile, condition of Theorem 5 does not satisfied, which means that the predator224 population also persist. However, one can see that the destabilization of the system has225 the potential to affect its persistence.226 227 M. Mohamed et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6097 15 of 19 0 20 40 60 80 100 0 1 2 3 4 5 Time P o p u la ti o n Prey Predator (a) Time series 0 1 2 3 4 5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 Prey P re d a to r (b) Phase plane Figure 4: Dynamical behaviour for the stochastic model (4)-(6) with initial values (u0 = 0.9, v0 = 0.8.), (δ = 0.4, ϵ = 1.2) and parameters values are given in (52). Figure 4 represents the stochastic model (4)-(6) with (δ = 0.4, ϵ = 1.2), in which the228 calculations shows the following:229 • ρu∗ k − αβu∗v∗ (1+Bu∗)2 − δ2 2 ≈ −0.0363 < 0.230 • λv∗ − ϵ2 2 ≈ 0− 0.6704 < 0.231 • ρ− δ2 2 ≈ 0.32 > 0.232 • γk ( ρ− δ2 2 ) − ρ ( δ + ϵ2 2 ) ≈ −0.288 < 0.233 From the above calculations, one can see that the conditions of the mean square stochas-234 tic stability Theorem 2 are not satisfied, so the system is unstable. Also, conditions of235 Theorems 3 and 5 are satisfied, which imply that the prey population u is weakly persist236 in average, while the predator population v go to extinct.237 238 0 20 40 60 80 100 0 1 2 3 4 5 Time P o p u la ti o n Prey Predator (a) Time series 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 Prey P re d a to r (b) Phase plane Figure 5: Dynamical behaviour for the stochastic model (4)-(6) with initial values (u0 = 0.9, v0 = 0.8.), (δ = 1.0, ϵ = 0.1) and parameters values are given in (52). M. Mohamed et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6097 16 of 19 Figure 5 represents the stochastic model (4)-(6) with (δ = 1.0, ϵ = 0.1), in which the239 calculations shows the following:240 • ρu∗ k − αβu∗v∗ (1+Bu∗)2 − δ2 2 ≈ −0.4563 < 0.241 • λv∗ − ϵ2 2 ≈ 00.0447 > 0.242 • ρ− δ2 2 ≈ −0.1 < 0.243 • γk ( ρ− δ2 2 ) − ρ ( δ + ϵ2 2 ) ≈ −0.723 < 0.244 Here, the calculations show that the conditions of Theorems 2 and 3 are not satisfied, while245 the condition of Theorem 5 is satisfied, which means that neither stability nor persistence246 of the stochastic system took place, but the system collapsed for both prey and predator247 populations. In this figure, high prey immigration results in species extinction, which248 contradicts the immigration argument. This phenomenon illustrates a case of destabiliza-249 tion within the ecological system due to the stochastic intensity of immigration, which250 prevents the system from self-regulating. In practice, this is referred to as the ”paradox251 of enrichment”.252 7. Conclusion253 In this paper, we investigated the effects of small random immigration in competitive254 predator- prey model with Holling type II. We obtained boundedness of the model’s solu-255 tion using both probabilistic and analytics inequalities. The conditions of stochastic mean256 square stability were derived (Theorem 2 and Corollary 1). A series of results on persis-257 tence and extinction for both prey and predator populations was obtained (Theorems 3 -258 5 and Corollary 2). All theoretical results were numerically verified in Figures 1-5. From259 this work, we conclude the following:260 • For fixed values of the parameters set (52), the stochastic stability of the model de-261 pends on the intensities of random immigration. When random immigration is small262 enough, the system becomes stochastically stable. However, as random immigration263 increases, the system is destabilized with high oscillations.264 • The persistence of the prey population u clearly affected by δ; where δ was small265 enough, that is, ρ− δ2 2 , the prey population u would persist with probability almost266 surely. As in condition (50) when ϵ was large enough, the predator population will267 extinct with probability almost surely. However, conditions (45) and (50)together268 emphasize that when δ was large enough (the immigration on the prey population269 was large enough), the whole system would extinct. This case is demonstrated in270 Figure 5.271 To the best of the authors’ knowledge, the first study to examine a competitive stochastic272 predator-prey model under the influence of small random immigration was [22], where273 M. Mohamed et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6097 17 of 19 the analysis focused on stochastic stability using a simpler model with a Holling type I274 functional response. In this paper, we present a comprehensive study of a more general275 and realistic competitive stochastic predator-prey model incorporating a Holling type II276 functional response, also subject to small random immigration. We investigate stochastic277 mean square stability and establish conditions for persistence and extinction of the species.278 In addition, our theoretical results are reinforced through numerical simulations, which279 illustrate and validate the analytical findings.280 Acknowledgements281 The authors extend the appreciation to the Deanship of Postgraduate Studies and282 Scientific Research at Majmaah University for funding this research work through the283 project number (ICR-2025-1724).284 References285 [1] Alfred James Lotka. Elements of physical biology. Williams & Wilkins, 1925.286 [2] Vito Volterra. Variazioni e fluttuazioni del numero d’individui in specie animali287 conviventi. Società anonima tipografica” Leonardo da Vinci”, 1926.288 [3] Patrick H Leslie. Some further notes on the use of matrices in population mathematics.289 Biometrika, 35(3/4):213–245, 1948.290 [4] Crawford S Holling. Some characteristics of simple types of predation and parasitism1.291 The canadian entomologist, 91(7):385–398, 1959.292 [5] Michael L Rosenzweig and Robert H MacArthur. Graphical representation and stabil-293 ity conditions of predator-prey interactions. The American Naturalist, 97(895):209–294 223, 1963.295 [6] Robert M May. Limit cycles in predator-prey communities. Science, 177(4052):900–296 902, 1972.297 [7] Herbert I Freedman. Deterministic mathematical models in population ecology. M.298 Dekker, 1980.299 [8] Yang Kuang. Basic properties of mathematical population models. J. Biomath,300 17(2):129–142, 2002.301 [9] Seda Igret Araz and Salah Boulaaras. Fractional modelling of gradual incorporation of302 infected prey into the predator-prey system with consideration of seasonality. Applied303 Mathematics in Science and Engineering, 33(1):2483728, 2025.304 [10] Zhihui Ma, Shufan Wang, Weide Li, and Zizhen Li. The effect of prey refuge in a305 patchy predator–prey system. Mathematical Biosciences, 243(1):126–130, 2013.306 [11] Sahabuddin Sarwardi, Prashanta Kumar Mandal, and Santanu Ray. Analysis of a307 competitive prey–predator system with a prey refuge. Biosystems, 110(3):133–148,308 2012.309 [12] Xiaoying Wang, Liana Zanette, and Xingfu Zou. Modelling the fear effect in predator–310 prey interactions. Journal of mathematical biology, 73(5):1179–1204, 2016.311 M. Mohamed et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6097 18 of 19 [13] Hang Deng, Fengde Chen, Zhenliang Zhu, and Zhong Li. Modelling the fear effect in312 predator–prey interactions. Advances in Difference Equations, 2019:1–17, 2019.313 [14] Takeru Tahara, Maica Krizna Areja Gavina, Takenori Kawano, Jerrold M Tubay,314 Jomar F Rabajante, Hiromu Ito, Satoru Morita, Genki Ichinose, Takuya Okabe,315 Tatsuya Togashi, et al. Asymptotic stability of a modified Lotka-Volterra model with316 small immigrations. Scientific reports, 8(1):7029, 2018.317 [15] Debasis Mukherjee. The effect of refuge and immigration in a predator–prey sys-318 tem in the presence of a competitor for the prey. Nonlinear Analysis: Real World319 Applications, 31:277–287, 2016.320 [16] Jawdat Alebraheem. Dynamics of a predator–prey model with the effect of oscillation321 of immigration of the prey. Diversity, 13(1):23, 2021.322 [17] Leonid Shaikhet and Andrei Korobeinikov. Persistence and Stochastic Extinction323 in a Lotka–Volterra Predator–Prey Stochastically Perturbed Model. Mathematics,324 12(10):1588, 2024.325 [18] Qiuyue Zhao and Xinglong Niu. Dynamics of a Stochastic Predator–Prey Model with326 Smith Growth Rate and Cooperative Defense. Mathematics, 12(12):1796, 2024.327 [19] Feng Rao and Yun Kang. Dynamics of a stochastic prey–predator system with328 prey refuge, predation fear and its carry-over effects. Chaos, Solitons & Fractals,329 175:113935, 2023.330 [20] Jaouad Danane, Karam Allali, Zakia Hammouch, and Kottakkaran Sooppy Nisar.331 Mathematical analysis and simulation of a stochastic COVID-19 Lévy jump model332 with isolation strategy. Results in Physics, 23:103994, 2021.333 [21] Mudhafar F Hama, Rando RQ Rasul, Zakia Hammouch, Kawa AH Rasul, and Jaouad334 Danane. Analysis of a stochastic SEIS epidemic model with the standard Brownian335 motion and Lévy jump. Results in Physics, 37:105477, 2022.336 [22] Jawdat Alebraheem, Mogtaba Mohammed, Ismail M Tayel, and Muhamad337 Hifzhudin Noor Aziz. Stochastic prey-predator model with small random immigra-338 tion. AIMS Mathematics, 9(6):14982–14996, 2024.339 [23] Hamlet Castillo-Alvino and Marcos Marvá. The competition model with Holling340 type II competitive response to interfering time. Journal of Biological Dynamics,341 14(1):222–244, 2020.342 [24] S Toaha et al. Stability analysis of prey predator model with holling ii functional343 response and threshold harvesting for the predator. In Journal of Physics: Conference344 Series., page 062025. IOP Publishing, 2019.345 [25] Jawdat Alebraheem. Relationship between the paradox of enrichment and the dy-346 namics of persistence and extinction in prey-predator systems. Symmetry, 10(10):532,347 2018.348 [26] Jawdat Alebraheem, Tabarek Qasim Ibrahim, Ghassan Ezzulddin Arif, Aws Asaad349 Hamdi, Omar Bazighifan, and Ali Hasan Ali. The stabilizing effect of small prey350 immigration on competitive predator-prey dynamics. Mathematical and Computer351 Modelling of Dynamical Systems, 30(1):605–625, 2024.352 [27] Doina Cioranescu, Patrizia Donato, and Marian P Roque. An introduction to sec-353 ond order partial differential equations: Classical and variational solutions. World354 M. Mohamed et al. / Eur. J. Pure Appl. Math, 18 (2) (2025), 6097 19 of 19 Scientific, 2018.355 [28] Valerii Nikolaevich Afanasiev, V Kolmanovskii, and Valerij Romanovič Nosov. Math-356 ematical theory of control systems design. Springer Science & Business Media, 2013.357 [29] Youlin Huang, Wanying Shi, Chunjin Wei, and Shuwen Zhang. A stochastic predator–358 prey model with Holling II increasing function in the predator. Journal of Biological359 Dynamics, 15(1):1–18, 2021.360 [30] Nirav Dalal, David Greenhalgh, and Xuerong Mao. A stochastic model for internal361 HIV dynamics. Journal of Mathematical Analysis and Applications, 341(2):1084–362 1101, 2008.363 Introduction Problem Statement Model Boundedness Stochastic mean square stability Persistence and extinction Persistence of the prey population Extinction of the prey population Extinction for the predator population Numerical Simulations and discussion Conclusion