Investigating Stochastic Dynamics of the Gilpin-Ayala Model in

Dispersed Polluted Environments

A. NAIT BRAHIM, B. HARCHAOUI, S. BOUTOUIL, M. EL IDRISSI, S. AZNAGUE,

A. SETTATI, A. LAHROUZ, M. EL JARROUDI

Department of Mathematics and Applications

Abdelmalek Essaadi University

Laboratory of Mathematics and Applications, FSTT, Abdelmalek Essaadi University, Tetouan, Morocco

MOROCCO

Abstract: This research delves into the analysis of a stochastic Gilpin-Ayala model operating within an anxious

environment, encompassing the phenomenon of diffusion between two distinct and specified geographical regions

that are the subjects of investigation. Initially, we rigorously formulate the essential criteria for ascertaining the

survival or extinction of the species. Furthermore, we furnish empirical substantiation for the presence of a stable

distribution. A significant milestone of our study involves the discernment and comprehensive delineation of

the pivotal determinants that intricately regulate extinction dynamics and persistence within the framework of

pollution parameters. This outcome underscores the pronounced impact of pollution on ecological dynamics

and affirms the necessity of incorporating pollution parameters into the purview of environmental investigations.

This revelation demonstrates that in the absence of pollution, the conventional criteria governing extinction and

persistence closely parallel those witnessed in unpolluted environments, thus validating the robustness of our

mathematical analysis. A series of numerical depictions are introduced to validate and provide empirical support

for the acquired results.

Key-Words: Stochastic Models, White Noise, Environmental Pollution, Extinction, Persistence, Stationary

Distribution, Lyapunov Function, Threshold

Received: January 13, 2023. Revised: July 14, 2023. Accepted: August 8, 2023. Published: September 13, 2023.

1 Introduction

Investigating ecological systems in the presence of

environmental contaminants has become a crucial and

highly regarded field of scientific research. There is

a growing concern regarding the impact of pollution

on biodiversity and the overall health of ecosystems.

Within this realm, it is particularly significant to com-

prehend the behavior of populations in environments

that are dispersed with pollutants. To address this,

the Gilpin-Ayala model has emerged as an invaluable

tool for studying the effects of pollutants on popula-

tion dynamics. This model allows researchers to gain

valuable insights into the interplay between species

interactions and environmental stressors. Population

ecology is a specialized branch of ecology that fo-

cuses on studying the intricate dynamics of species

populations and their interactions with the environ-

ment. However, due to the accelerated growth of in-

dustries and agriculture, a wide range of toxic sub-

stances are being discharged into the atmosphere, pos-

ing significant risks to the survival of exposed or-

ganisms. The proliferation of species in the natural

realm is intricately interconnected with the influence

of environmental noise (refer to [1], [2], [3], [4], [5],

[6], [7], [8], [9], for further details). Consequently,

several esteemed authors have delved into stochas-

tic population models to understand the dynamics

within polluted environments (see, for instance, [10],

[11]). Notably, [12] conducted an in-depth investiga-

tion into a stochastic model that specifically examines

the effects of toxins on a single-species system.











dx =x(t) [r0−k0x(t)−l1c0(t)] dt

+σ1x(t)dB1(t),

dc0= [−(g+m)c0(t) + kce(t)] dt,

dce= [−hce(t) + u(t)] dt,

(1)

where x(t)represents the population size at time t, the

parameter r0>0corresponds to the inherent growth

rate of the population in the absence of any toxicant.

Furthermore, l1>0indicates the response rate of

the population to the pollutant present in the organ-

ism. The variables c0and cedenote the concentrations

of noxious compounds within the organism and the

surrounding environment, respectively. Additionally,

B1denotes a conventional standard Brownian mo-

tion, delineated within the confines of a comprehen-

sive probability space denoted as (Ω,F,P). The pa-

rameter σ1signifies the magnitude of the white noise

intensity, while k > 0signifies the rate at which the

organism absorbs the toxicant from its surrounding

environment. The parameters g > 0and m > 0rep-

resent the egestion and depuration rates of the toxicant

WSEAS TRANSACTIONS on MATHEMATICS

DOI: 10.37394/23206.2023.22.67

A. Nait Brahim, B. Harchaoui, S. Boutouil,

M. El Idrissi, S. Aznague, A. Settati,

A. Lahrouz, M. El Jarroudi

E-ISSN: 2224-2880

607

Volume 22, 2023

within the organism, respectively. Similarly, the pa-

rameter h > 0represents the rate at which the toxicant

undergoes environmental loss. Furthermore, u(t)is a

non-negative, bounded, and continuous function de-

fined over the interval [0,+∞), representing the ex-

ogenous rate at which toxicants are introduced into

the environment. The authors [12] have effectively

ascertained the threshold governing species persis-

tence or extinction and concurrently established req-

uisite conditions for the stochastic permanence of the

population. Moreover, in a related study, [13], the

authors examined a scenario where both the param-

eters l1and k0in the model (1) are also subject to

random fluctuations. By applying the Gilpin-Ayala

model, they derived the following stochastic single-

species model, taking into account the effects of the

toxicant.











dx =x(t)[r0−k0xθ1(t)−l1c0(t)]dt

+σ1x(t)dB1(t) + σ2c0(t)x(t)dB1(t)

+σ1x1+θ1(t)dB3(t),

dc0= [ −(g+m)c0(t) + kce(t)]dt,

dce= [ −hce(t) + u(t)]dt.

(2)

However, species dispersal is a widely recognized oc-

currence in the natural world. It frequently transpires

within patches of ecological environments, mainly

due to the influence of human activities and indus-

tries on the environment. Various factors have par-

titioned reproductive and population-centric domains

and other habitats into discrete patches. These factors

encompass the spatial distribution of industrial facil-

ities and the pollution of the air, soil, and water bod-

ies. This issue has been extensively studied and doc-

umented (refer to [14], [15], [16], [17]). The present

study aims to examine the implications of the disper-

sal phenomenon. To achieve this, we have developed

a stochastic diffusion system comprising two patches

affected by a toxicant.











dx1=x1r1−k1xθ1

1−l1c0(t)

+ε12(x2−x1)dt +

i=1 α1ix1

+β1ix1+θ1

1+γ1ix1c0(t)dBi,

dx2=x2r2−k2xθ2

2−l2c0(t)

+ε21(x1−x2)dt +

i=1 α2ix2

+β2ix1+θ2

2+γ2ix1c0(t)dBi,

dc0= [−(g+m)c0(t) + kce(t)] dt,

dce= [−hce(t) + u(t)] dt.

(3)

Within this framework, xisignifies the population

density for a particular species residing in the ith

patch. Similarly, the riand kidenote the growth rate

and self-competition coefficient that regulate the pop-

ulation dynamics within those above the ith patch.

Furthermore, we assign the symbol εi,j >0to rep-

resent a diffusion coefficient. This represents a non-

negative parameter that characterizes species migra-

tion from the jth patch to the ith patch (with the con-

dition i=j). We assume that the net migration of

individuals from the jth patch to the ith patch is di-

rectly proportional to the disparity in population den-

sities (xi−xj) between these patches. This assump-

tion aligns with the commonly accepted framework

in similar research endeavors (refer to [18], [19], for

more details). The vectors αi= (α1i, α2i, ..., αni)

and βi= (β1i, β2i, ..., βni)represent the magni-

tudes of white noise signals at positions riand ki,

respectively. In order to capture the correlation be-

tween the noises at riand ki, we utilize a vector

B(t) = (B1(t), B2(t), ..., Bn(t))T, which denotes an

n-dimensional Brownian motion. In this study, we

consider a comprehensive probability space denoted

as (Ω,F,{Ft}t≥0,P), where Ωrepresents the sam-

ple space, Fdenotes the sigma-algebra of events, and

{Ft}t≥0is a filtration satisfying the standard condi-

tions. Throughout our analysis, we use the custom-

ary inner product ⟨., .⟩and the Euclidean norm |.|de-

fined on Rn. The paper is structured as follows: Sec-

tion 2 addresses the issue of determining the existence

and uniqueness of a positive solution for the system

described in equation (1). In Section 3, a rigorous

survival analysis is conducted to establish sufficient

conditions for various ecological outcomes, encom-

passing scenarios of extinction, non-persistence in the

mean, weak persistence, and stochastic permanence

of the species. Section 4 further contributes to the

field by demonstrating the existence of a stationary

distribution, building upon and improving the previ-

ously proposed sufficient condition by [19]. Finally,

the main findings of the study are illustrated through

numerical simulations, offering visual insights into

the observed phenomena.

2 Existence and Uniqueness

With reference to [12], we have the following lemma.

Lemma 2.1 If

lim sup

t→∞

u(t)≤hand 0< k ≤g+m,

then, for each t≥0,

0≤ce(t)<1,0≤c0(t)<1.

Therefore, we assume that

0< k ≤g+mand lim sup

t→∞

u(t)≤h.

WSEAS TRANSACTIONS on MATHEMATICS

DOI: 10.37394/23206.2023.22.67

A. Nait Brahim, B. Harchaoui, S. Boutouil,

M. El Idrissi, S. Aznague, A. Settati,

A. Lahrouz, M. El Jarroudi

E-ISSN: 2224-2880

608

Volume 22, 2023

It is imperative to acknowledge that the final two

equations in a model (1) exhibit linearity concern-

ing c0(t)and ce(t). Consequently, deriving their

explicit solutions becomes straightforward. Subse-

quently, our focus narrows down to exclusively ad-

dressing the initial two equations within a model (1),

namely:











dx1=x1r1−l1c0(t)−k1xθ1

1

+ε12(x2−x1)dt +

i=1 α1ix1

+β1ix1+θ1

1+γ1ix1c0(t)dBi,

dx2=x2r2−l2c0(t)−k2xθ2

2

+ε21(x1−x2)dt +

i=1 α2ix2

+β2ix1+θ2

2+γ2ix1c0(t)dBi.

(4)

The population densities, denoted as x1and x2, pos-

sess biological significance, requiring them to be non-

negative. To address this requirement, we will ana-

lyze system (4) within a specific region

+={(x1, x2)|xi>0, i = 1,2}.

We will now demonstrate that the set R2

+is a positive

invariant set.

Theorem 2.2 For each (x1(0), x2(0)) ∈R2

+, there

exists a unique solution (x1(t), x2(t)) to system (4)

for t≥0. Furthermore, the solution will remain

within R2

+with a probability 1.

Proof 1 Since all the coefficients in system (4) ex-

hibit local Lipschitz continuity. Therefore, for any ini-

tial value (x1(0), x2(0)) ∈R2

+, there exists a unique

local solution (x1(t), x2(t)) for t∈[0, τe[, where

τerepresents the explosion time (refer to [20], [21],

[22], [23], for more details). To prove that this solu-

tion is global, it is necessary to prove that τe=∞.

Let k0≥0be such that

xi(0) ∈1

, k0, i = 1,2.

For each integer k≥k0, we define stopping times as

follows

τk= inf t∈[0, τe]; xi(t)/∈1

, k, i = 1,2.

Set τ∞= lim

k→+∞τk, thus τ∞≤τea.s..

Assume that τ∞finite, then there are two constants

T > 0and ε∈(0,1), where P(τ∞≤T)>2ε.

Then, there exists k1≥k0, where P(τk≤T)≥ε.

There is an integer k≥k1, denote Ωk={τk≤T},

then

P(Ωk)≥ε. (5)

Define the function V∈ C2(R2

+;R+)such that

V(x1, x2) = [2√x1−ln(x1)]

+ [2√x2−ln(x2)] .(6)

It is important to note that for each ω∈Ωk, there ex-

ists an isuch that xi(τk, ω)takes the value of either k

or 1

k. Moreover, the function V(x1(τk, ω), x2(τk, ω))

is guaranteed to be greater than or equal to either

h√k−1−0.5 ln(k)ior 1

√k−1−0.5 ln 1

k.

Hence, we can deduce

V(x1(τk, ω), x2(τk, ω)) ≥ˇ

V(k),(7)

where

V(k) = h√k−1−0.5 ln(k)i

∧1

√k−1−0.5 ln 1

k.

Applying Itô formula, we obtain

dV =x−0.5

1−x−1

1x1r1−l1c0(t)−k1xθ1

1

+ε12 (x2−x1)dt +

i=1 α1ix1+β1ix1+θ1

+γ1ix1c0(t)dBi+x−0.5

2−x−1

2

×x2(r2−l2c0(t)−k2xθ2

+ε21(x1−x2)dt +

i=1 α2ix2+β2ix1+θ2

+γ2ix1c0(t)dBi+−0.25x−1.5

1+ 0.5x−2

1

i=1

1α1i+β1ixθ1

1+γ1ic0(t)2

+−0.25x−1.5

2+ 0.5x−2

2

i=1

2α2i

+β2ixθ2

2+γ2ic0(t)2

dt.

WSEAS TRANSACTIONS on MATHEMATICS

DOI: 10.37394/23206.2023.22.67

A. Nait Brahim, B. Harchaoui, S. Boutouil,

M. El Idrissi, S. Aznague, A. Settati,

A. Lahrouz, M. El Jarroudi

E-ISSN: 2224-2880

609

Volume 22, 2023

For i= 1,2, there exists Ni>0, where

x−0.5

i−x−1

i< Ni,for t > 0.

Letting N= max{N1, N2}, we get

dV ≤N−k1x1+θ1

1+ (r1−l1c0(t) + ε21

−ε12)x1dt +−0.25x0.5

1+ 0.5

i=1 α1i

+β1ixθ1

1+γ1ic02

dt +N−k2x1+θ2

+(r2−l2c0(t) + ε12 −ε21)x2dt

+−0.25x0.5

2+ 0.5

i=1 α2i+β2ixθ2

+γ2ic02

dt + (x0.5

1−1)

i=1 α1i+β1ixθ1

+γ1ic0(t)dBi+x0.5

2−1

i=1 α2i

+β2ixθ2

2+γ2ic0(t)dBi.(8)

Denote

g(x1) = N−k1x1+θ1

1+ (r1−l1c0(t) + ε21

−ε12)x1dt +−0.25x0.5

1+ 0.5

i=1 α1i+β1ixθ1

1+γ1ic0(t)2

dt,

and

h(x2) = N−k2x1+θ2

2+ (r2−l2c0(t) + ε12

−ε21)x2dt + ( −0.25x0.5

2+ 0.5)

i=1 α2i+β2ixθ2

2+γ2ic0(t)2dt.

One can easily verify that

lim

x1→+∞g(x1) = lim

x2→+∞h(x2) = −∞,

and then there is M > 0such that

g(x1) + h(x2)< M. (9)

From (8) and (9), we obtain

dV ≤Mdt +x0.5

1−1

i=1 α1i+β1ixθ1

+γ1ic0(t)dBi(t) + x0.5

2−1

i=1 α2i+β2ixθ2

2+γ2ic0(t)dBi(t).

By integrating both sides from 0to τk∧T, one has

Zτk∧T

dV ≤Zτk∧T

Mdt +Zτk∧T

0x0.5

1−1

i=1 α1i+β1ixθ1

1+γ1ic0(t)dBi

+Zτk∧T

0x0.5

2−1

i=1 α2i

+β2ixθ2

2+γ2ic0(t)dBi(t).

Consequently

E[V(x1(τk∧T), x2(τk∧T))]

≤V(x1(0), x2(0)) + ME[(τk∧T)] ,

≤V(x1(0), x2(0)) + MT,

which yields

V(x1(0), x2(0)) + MT

≥E[1ΩkV(x1(τk, ω), x2(τk, ω))] .(10)

By (5) and (10), we obtain

V(x1(0), x2(0)) + MT

≥εV (x1(τk, ω), x2(τk, ω)) .(11)

From (7) and (11), one obtains

V(x1(0), x2(0)) + MT ≥εˇ

V(k).

Letting k−→ ∞, we get the following contradiction

V(x1(0), x2(0)) + MT =∞.

Therefore, we need to have τ∞=∞a.s..

3 Extinction

To comprehensively explore the intricate topic of

species extinction, [24], utilizing a crucial lemma

termed the exponential martingale inequality is

paramount (refer to [25], Theorem 7.4, page 44).

WSEAS TRANSACTIONS on MATHEMATICS

DOI: 10.37394/23206.2023.22.67

A. Nait Brahim, B. Harchaoui, S. Boutouil,

M. El Idrissi, S. Aznague, A. Settati,

A. Lahrouz, M. El Jarroudi

E-ISSN: 2224-2880

610

Volume 22, 2023

Lemma 3.1 Consider a local martingale (Mt)t≥0,

which vanishes at time zero. Thus, for each positive

constants a,b, and c, we have

Psup

0≤t≤aMt−b

2[Mt, Mt]> c≤exp (−bc),

where [Mt, Mt]is the quadratic variation of Mt.

Theorem 3.2 For all (x1(0), x2(0)) ∈R2

+, the solu-

tion of the SDE (4) obeys

lim sup

t→∞

tlog x1(t)

ε12

+x2(t)

ε21 ≤M−m2

2a.s.,

where

M= max r1−l1inf

t≥0{c0(t)};r2−l2inf

t≥0{c0(t)},

(12)

and

m= min

1≤i≤nα1i+γ1iinf

t≥0{c0(t)};α2i

+γ2iinf

t≥0{c0(t)}.

Moreover, if M−m2

2<0, then the species in (4)

is extinct.

Proof 2 We define

X(t) = x1(t)

ε12

+x2(t)

ε21

Using the Itô formula, we obtain

dlog(X) = 1

X(t)x1(t)

ε12

(r1−l1c0(t))

−k1xθ1

1(t) + x2(t)

ε21

(r2−l2c0(t))

−k2xθ2

2(t)dt −1

2 (X(t))2

i=1 x1(t)

ε12 α1i+β1ixθ1

1(t)

+γ1ic0(t)+x2(t)

ε21 α2i+β2ixθ2

2(t)

+γ2ic0(t)2

dt +1

X(t)

i=1 x1(t)

ε12

×α1i+β1ixθ1

1(t) + γ1ic0(t)

+x2(t)

ε21 α2i+β2ixθ2

2(t)

+γ2ic0(t)dBi.(13)

By integrating, we obtain

log (X(t)) = Zt

X(s)x1(s)

ε12

(r1−l1c0(s))

−k1xθ1

1(s) + x2(s)

ε21

(r2−l2c0(s))

−k2xθ2

2(s)ds −Zt

2 (X(s))2

i=1 x1(s)

ε12 α1i+β1ixθ1

1(s)(14)

+γ1ic0(s)+x2(s)

ε21 α2i+β2ixθ2

2(s)

+γ2ic0(s)2

ds +Mt+ log(X(0)),

where

Mt=Zt

X(s)

i=1 x1(s)

ε12 α1i+β1ixθ1

1(s)

+γ1ic0(s)+x2(s)

ε21 α2i+β2ixθ2

2(s)

+γ2ic0(s)dBi(s),

is a real-valued continuous local martingale vanish-

ing at t= 0 with the quadratic variation

[Mt, Mt] =Zt

X−2(s)

i=1 x1(s)

ε12 α1i+β1ixθ1

1(s)

+γ1ic0(s)+x2(s)

ε21 α2i+β2ixθ2

2(s)

+γ2ic0(s)2

ds. (15)

By Lemma 3.1, for any integer k≥1and ϵsufficiently

small, we get

P"sup

0≤t≤kMt−ϵ

2[Mt, Mt]>2

ϵlog(k)#≤1

k2.

By the Borel-Cantelli Lemma, there is Ω1⊂Ωwith

P(Ω1) = 1, and for each ω∈Ω1, there is k1(ω)such

that

Mt≤2

ϵlog(k) + ϵ

2[Mt, Mt],(16)

where, 0≤t≤kand k≥k1(ω).

WSEAS TRANSACTIONS on MATHEMATICS

DOI: 10.37394/23206.2023.22.67

A. Nait Brahim, B. Harchaoui, S. Boutouil,

M. El Idrissi, S. Aznague, A. Settati,

A. Lahrouz, M. El Jarroudi

E-ISSN: 2224-2880

611

Volume 22, 2023

Hence, from (14) and (16), we get for ω∈Ω1,

k≥k1(ω)and 0≤t≤k

log (X(t))

⩽Zt

X(s)x1(s)

ε12 r1−l1c0(s)

−k1xθ1

1(s)+x2(s)

ε21 r2−l2c0(s)

−k2xθ2

2(s)ds −Zt

1−ϵ

2 (X(s))2

i=1 x1(s)

ε12

(α1i+β1ixθ1

1(s) + γ1ic0(s))

+x2(s)

ε21 α2i+β2ixθ2

2(s) + γ2ic0(s)2

+ log (X(0)) + 2

ϵlog(k),(17)

which implies

log (X(t))

≤Zt

X(s)x1(s)

ε12

(r1−l1c0(s))

+x2(s)

ε21

(r2−l2c0(s)) ds −Zt

1−ϵ

2x1(s)

ε12

+x2(s)

ε21

t−2n

i=1 x1(s)

ε12

(α1i+γ1ic0(s))

+x2(s)

ε21

(α2i+γ2ic0(s))2

ds + log (X(0))

ϵlog(k).(18)

Consequently, it is evident from (18) that

log (X(t)) ≤M−1−ϵ

2m2t+ log (X(0))

ϵlog(k).(19)

Take into account ω∈Ω1and consider a value of t

that is sufficiently large such that the maximum inte-

ger smaller than tfulfills the condition [t]≥k1(ω).

Based on equation (19), we get

tlog (X(t)) ≤M−1−ϵ

2m2+1

[t]log (X(0))

ϵ(2 log ([t] + 1) ),(20)

which gives

lim sup

t→∞

tlog (X(t)) ≤M−1−ϵ

2m2.

Letting ϵ−→ 0, one obtains

lim sup

t→∞

tlog (X(t)) ≤M−1

2m2.

This proves the theorem.

Theorem 3.3 If β1i=β2i= 0, then

(i) lim inf

t→∞

tZt

xθ1

1(s)ds ≥1

k1r1−l1−ε12

−1

2∥α1+γ1∥2,

(ii) lim inf

t→∞

tZt

xθ2

2(s)ds ≥1

k2r2−l2−ε21

−1

2∥α2+γ2∥2.

Proof 3 (i)From (2), β1i=β2i= 0 and Itô formula,

we get

dlog(x1(t)) ≥r1−l1−k1xθ1

1(t)+ε12

x1(t)

×(x2(t)−x1(t)) −1

2∥α1+γ1∥2dt

i=1

(α1i+γ1ic0(t)) dBi,(21)

which implies

dlog(x1(t))

≥r1−l1−ε12 −1

2∥α1+γ1∥2(22)

−k1xθ1

1(t)dt +

i=1

(α1i+γ1ic0(t)) dB.

Integrating, we get

log (x1(t)) ≥r1−l1−ε12 −1

2∥α1+γ1∥2t

−k1Zt

xθ1

1(s)ds +Zt

i=1

(α1i

+γ1ic0(t))dBi+ log(x1(0)).

WSEAS TRANSACTIONS on MATHEMATICS

DOI: 10.37394/23206.2023.22.67

A. Nait Brahim, B. Harchaoui, S. Boutouil,

M. El Idrissi, S. Aznague, A. Settati,

A. Lahrouz, M. El Jarroudi

E-ISSN: 2224-2880

612

Volume 22, 2023

Hence

log xθ1

1≤θ1r1−l1−ε12 −1

2∥α1+γ1∥2t

−θ1k1Zt

xθ1

1(s)ds

+θ1Zt

i=1

(α1i+γ1ic0(t))dBi

+ log(xθ1

1(0)).(23)

Thus

xθ1

1(t) exp θ1k1Zt

xθ1

1(s)ds

≥exp θ1r1−l1−ε12 −1

2∥α1+γ1∥2t

+θ1Nt+ log xθ1

1(0),(24)

where

Nt=Zt

i=1

(α1i+γ1ic0(t)) dBi,

represents a continuous martingale with real-valued

properties characterized by its quadratic variation.

[Nt, Nt] = Zt

i=1

(α1i+γ1ic0(t))2ds.

According to the law of large numbers for martingales

(refer to [25], for further information), one obtains

lim

t→∞

tθ1r1−l1−ε12 −1

2∥α1+γ1∥2t

+θ1Nt+ log xθ1

1(0)(25)

=θ1r1−l1−ε12 −1

2∥α1+γ1∥2a.s..

From (25), one can easily show that there exists Ω′

1⊂

Ωsuch that P(Ω′

1)=1and for each ω∈Ω′

1and

ϵ > 0sufficiently small, there is T(ω, ϵ)such that for

any t≥T, we have



tθ1r1−l1−ε12 −1

2∥α1+γ1∥2t+θ1Nt

+ log xθ1

1(0)−θ1r1−l1−ε12

−1

2∥α1+γ1∥2≤ϵ.

θ1r1−l1−ε12 −1

2∥α1+γ1∥2t+θ1Nt

+ log xθ1

1(0)

≥θ1r1−l1−ε12 −1

2∥α1+γ1∥2−ϵt,

which gives with (24) that for any t≥T

xθ1

1(t) exp θ1k1Zt

xθ1

1(s)ds(26)

≥exp θ1r1−l1−ε12 −1

2∥α1+γ1∥2−ϵt,

xθ1

1(t) exp θ1k1Zt

xθ1

1(s)ds(27)

θ1k1

dt exp θ1k1Zt

xθ1

1(s)ds.

By (26) and (27), we have

θ1k1

dt exp θ1k1Zt

xθ1

1(s)ds (28)

≥exp θ1r1−l1−ε12 −1

2∥α1+γ1∥2−ϵt.

Hence, by integrating (28) from Tto tyields that

θ1k1exp θ1k1Zt

xθ1

1(s)ds

−exp θ1k1ZT

xθ1

1(s)ds

≥1

θ1r1−l1−ε12 −1

2∥α1+γ1∥2−ϵ−1

×exp θ1r1−l1−ε12 −1

2∥α1+γ1∥2

−ϵt−exp θ1r1−l1−ε12

−1

2∥α1+γ1∥2−ϵT.

Then, we have

tZt

xθ1

1(s)ds ≥1

θ1k1tΛ(t),(29)

WSEAS TRANSACTIONS on MATHEMATICS

DOI: 10.37394/23206.2023.22.67

A. Nait Brahim, B. Harchaoui, S. Boutouil,

M. El Idrissi, S. Aznague, A. Settati,

A. Lahrouz, M. El Jarroudi

E-ISSN: 2224-2880

613

Volume 22, 2023

where

Λ(t)=log exp θ1k1ZT

xθ1

1(s)ds

+k1r1−l1−ε12 −1

2∥α1+γ1∥2−ϵ−1

×exp θ1r1−l1−ε12 −1

2∥α1+γ1∥2

−ϵt−exp θ1r1−l1−ε12

−1

2∥α1+γ1∥2−ϵT.

r1−l1−ε12 −1

2∥α1+γ1∥2−ϵ≤0,

then

lim

t→∞

tΛ(t) = 0.

Hence the assertion (i)holds trivially. Else if

r1−l1−ε12 −1

2∥α1+γ1∥2−ϵ>0,

we easily have

lim

t→∞

Λ(t)

t=θ1r1−l1−ε12 −1

2∥α1+γ1∥2−ϵ.

Using (29) and (30), we get

lim inf

t→∞

tZt

xθ1

1(s)ds ≥1

k1r1−l1−ε12

−1

2∥α1+γ1∥2−ϵ.

By letting ϵ−→ 0, we get the required estimation (i).

Similarly, we get the assertion (ii).

Based on Theorem 3.3, we can deduce the subsequent

corollary

Corollary 3.1 Suppose that β1i=β2i= 0, thus

(i)If r1−l1−ε12 −1

2∥α1+γ1∥2>0, then xθ1

is strongly persistent in mean,

(ii)If r2−l2−ε21 −1

2∥α1+γ2∥2>0, then xθ2

is strongly persistent in mean.

4 Stationary distribution

In this section, our investigation determines whether

a solution to the SDE represented by equation (4) ex-

hibits an asymptotically invariant distribution, indi-

cating stability in a stochastic context. To shed light

on this matter, we turn to the insightful analysis con-

ducted by [26], whose theorem provides a compre-

hensive answer to this question. Let us consider a ho-

mogeneous Markov process denoted as X(t), which

is characterized by the following stochastic differen-

tial equation

dX(t) = b(X)dt +

r=1

σr(X)dBr(t),(30)

X(t) = x1(t)

x2(t), b(X) = b1(x1, x2)

b2(x1, x2),

b(X) = 





x1(t)r1−l1c0(t)−k1xθ1

1(t)

+ε12 (x2(t)−x1(t))

x2(t)r2−l2c0(t)−k2xθ2

2(t)

+ε21 (x1(t)−x2(t))







σr(X) = σ1

r(x)

σ2

r(x),

σr(X) = 





α1rx1(t) + β1rx1+θ1

1(t)

+γ1rx1c0(t)

α2rx2(t) + β2rx1+θ2

2(t)

+γ2rx2c0(t)







Let V(x)∈ C2(R2)be a function that is twice con-

tinuously differentiable. The differential operator L

mentioned in equation (30) can be defined as follows

LV(x) = ∇V(x)b(x) + 1

2T r A(x)∇2V(x),

where the gradient of the function V(x), denoted

as ∇V(x), and the hessian of V(x), represented as

∇2V(x), play crucial roles in this context. Addition-

ally, the diffusion matrix Ais defined to be associated

with these quantities.

A(x) = (aij (x))1⩽i,j⩽2, aij =

r=1

σi

r(x)σj

r(x).

So, the diffusion matrix associated with the system

(30) can be expressed as follows

A(x1, x2) = ∥Y∥2< Y, Z >

< Y, Z > ∥Z∥2,(31)

WSEAS TRANSACTIONS on MATHEMATICS

DOI: 10.37394/23206.2023.22.67

A. Nait Brahim, B. Harchaoui, S. Boutouil,

M. El Idrissi, S. Aznague, A. Settati,

A. Lahrouz, M. El Jarroudi

E-ISSN: 2224-2880

614

Volume 22, 2023

where

Y=hα1rx1(t) + β1rx1+θ1

1(t) + γ1rx1c0(t)i1≤i≤n,

Z=hα2rx2(t) + β2rx1+θ2

2(t) + γ2rx2c0(t)i1≤i≤n.

We can readily ascertain the positive definiteness of

matrix Aby verifying the strict inequality of the

Cauchy-Schwarz inequality, expressed as

|< Y, Z > |≤∥Y∥.∥Z∥.

In cases where a system (30) lacks equilibrium states,

exploring the potential existence of an asymptotically

invariant distribution remains feasible. The authors

[26] propose that it is sufficient to establish the pos-

itive recurrence of system solutions (30) concerning

a bounded open set to pursue this investigation. We

consider the R2-valued process X(t, x0)is recurrent

for the bounded set U⊂R2if P(τx0<∞) = 1, for

each x0∈Uc. Here, τx0represents the stopping time

of Ufor the process X(t, x0), defined as

τx0= inf{t > 0, X(t, x0)∈U}.

The process X(t, x0)is said to be positive recurrent

for the set Uif it satisfies two conditions: first, it is

recurrent for U, meaning that it revisits the set Uin-

finitely often, and second, for any x0∈ U,E(τx0)<

∞. A theorem exists that provides a criterion for pos-

itive recurrence based on the Lyapunov function (see,

e.g., [26], and the references cited therein).

Theorem 4.1 The system represented by (2) is said to

be positively recurrent if there exists a bounded open

subset D⊂R2with a smooth boundary that satisfies

the following conditions

(i)There is some κ∈(0,1] such that for all

(x1, x2)∈D,ξ∈R2

κ∥ξ∥2≤ξTA(x1, x2)ξ≤κ−1∥ξ∥2,

(ii)There exists V∈ C(Dc;R+)a nonnegative func-

tion that is twice continuously differentiable and for

some ϱ > 0

LV(x1, x2)≤ −ϱfor all (x1, x2)∈Dc.

Furthermore, the system (2) possesses a distinctive

ergodic stationary distribution denoted as π, and the

solution (x1(t), x2(t)) exhibits uniqueness concern-

ing this distribution. If the function fis integrable

with respect to the measure π, then

Plim

t→∞

tZt

f(X(s))ds =ZR2

f(x)π(dx)= 1.

The subsequent theorem provides a condition suffi-

cient for a stationary distribution in our model (2).

Theorem 4.2 Consider the stochastic system (2) with

an initial condition in R2

+, where n≥4. Let us as-

sume that α1,α2,β1,β2,γ1, and γ2are linearly in-

dependent and satisfy the following conditions

r1−l1−ε12 −1

2∥α1+γ1∥2>0,(32)

and

r2−l2−ε21 −1

2∥α2+γ2∥2>0.

Then, the solution (x1(t), x2(t)) of the SDE (2) admits

a unique stationary distribution and is ergodic.

Proof 4 Given a bounded open subset as described

below

D=1

µ, µ×1

µ, µ⊂R2

+,(33)

where µrepresents a sufficiently large number. So,

D⊂R2

(i)If α1, α2, β1, and β2are linearly independent, it

follows that Yand Zare linearly independent. Con-

sequently, based on the abovementioned observation,

the matrix Ahas the potential to be positive and def-

inite. This leads us to the following expression:

λmin =λmin (A(x1, x2)) >0,(34)

and

λmax =λmax (A(x1, x2)) >0,

where the quantities λmin (A(x1, x2)) and

λmax (A(x1, x2)) denote the smallest and largest

eigenvalues of the matrix A(x1, x2), respectively.

In addition, we can establish the following for all

ξ∈R2

λmin∥ξ∥2⩽ξTA(x1, x2)ξ⩽λmax∥ξ∥2.(35)

It is evident that the functions λmin (A(., .)) and

λmax (A(., .)) are continuous with respect to the vari-

ables (x1, x2). Consequently, using equation (34), we

can deduce that

λ1= min

(x1,x2)∈D

λmin (A(x1, x2)) >0,(36)

and

λ2= max

(x1,x2)∈D

λmax (A(x1, x2)) >0.

WSEAS TRANSACTIONS on MATHEMATICS

DOI: 10.37394/23206.2023.22.67

A. Nait Brahim, B. Harchaoui, S. Boutouil,

M. El Idrissi, S. Aznague, A. Settati,

A. Lahrouz, M. El Jarroudi

E-ISSN: 2224-2880

615

Volume 22, 2023

Furthermore, by (35), we get, for all ξ∈R2

κ∥ξ∥2≤ξTA(x1, x2)ξ≤1

κ∥ξ∥2,

where κ= min{λ1, λ2,1}. Consequently, we can

confirm that the condition (i)specified in Theorem

4.1 holds for SDE (2).

(ii)Consider the following positive functions

ψ1(x1) = 1

2log2(x1), ψ2(x2) = 1

2log2(x2),

and

ψ3(x1, x2) = ε21x1+ε12x2,

and

ψ(x1, x2) = ψ1(x1) + ψ2(x2) + ψ3(x1, x2).

Using Itô formula on the function ψ1, we get

Lψ1(x1) = log(x1)r1−l1c0−k1xθ1

+ε12 x2

x1−1+1

2(1 −log(x1))

r=1 α1r+β1rxθ1

1+γ1rc02.

Using log(x1)≤x1and rearranging yields

Lψ1(x1)≤r1−l1c0−ε12 −1

2∥α1+γ1c0∥2

×log(x1) + ε12x2−k1xθ1

1log(x1)

2∥α1+γ1c0∥2(37)

22<(α1+γ1c0), β1>

+∥β1∥2xθ1

1xθ1

1(1 −log(x1)).

Similarly, we have

Lψ2(x2)≤r2−l2c0−ε21 −1

2∥α2+γ2c0∥2

×log(x2) + ε21x1−k2xθ2

2log(x2)

2∥α2+γ2c0∥2

22<(α2+γ2c0), β2>(38)

+∥β2∥2xθ2

2xθ2

2(1 −log(x2)) ,

and

Lψ3(x1, x2) = ε21 r1x1−l1c0x1−k1x1+θ1

1

+ε12 r2x2−l2c0x2−k2x1+θ2

2.

(39)

From (37), (38) and (39), we have

Lψ(x1, x2)≤ϕ1(x1) + ϕ2(x2),(40)

where

ϕ1(x1) = r1−l1c0−ε12 −1

2∥α1+γ1c0∥2

×log(x1) + 1

22<(α1+γ1c0), β1>

+∥β1∥2xθ1

1xθ1

1(1 −log(x1))

−k1xθ1

1log(x1)+(ε21r1+ε21

−ε21l2c0)x1−k1ε21x1+θ1

2∥α1+γ1c0∥2,

and

ϕ2(x2) = r2−l2c0−ε21 −1

2∥α2+γ2c0∥2

×log(x2) + 1

22<(α2+γ2c0), β2>

+∥β2∥2xθ2

2xθ2

2(1 −log(x2))

−k2xθ2

2log(x2)+(ε12r2+ε12

−ε12l1c0)x2−k2ε12x1+θ2

2∥α2+γ2c0∥2.

Hence, if x1∼ ∞, then

ϕ1(x1)≃1

2∥β1∥2x2θ1

1log(x1)−k1ε21x1+θ1

and if x1∼0, then

ϕ1(x1)≃r1−l1−ε12 −1

2∥α1+γ1∥2log(x1),

and if x2∼ ∞, then

ϕ2(x2)≃ −1

2∥β2∥2x2θ2

2log(x2)−k2ε12x1+θ2

and if x2−→ 0, then

ϕ2(x2)≃r2−l2−ε21 −1

2∥α1+γ2∥2log(x2).

WSEAS TRANSACTIONS on MATHEMATICS

DOI: 10.37394/23206.2023.22.67

A. Nait Brahim, B. Harchaoui, S. Boutouil,

M. El Idrissi, S. Aznague, A. Settati,

A. Lahrouz, M. El Jarroudi

E-ISSN: 2224-2880

616

Volume 22, 2023

Since

r1−l1−ε12 −1

2∥α1+γ1∥2>0,

and

r2−l2−ε21 −1

2∥α1+γ2∥2>0,

then

lim

x1→0ϕ1(x1) = lim

x1→+∞ϕ1(x1) = −∞,

and

lim

x2→0ϕ2(x2) = lim

x2→+∞ϕ2(x2) = −∞.

This gives together with (33) and (40) that for a suffi-

ciently large µ,

Lψ(x1, x2)≤ −1for each (x1, x2)∈Dc.

As a result, the stochastic system described by equa-

tion (2) possesses an invariant distribution, which is

characterized by a density of zero in R2

5 Computer simulations

To demonstrate our findings, we will utilize the

widely-known Euler scheme (refer to, [27], for fur-

ther details). We shall analyze the discretized system

presented below











x1(k+ 1) = x1(k) + x1(k)r1−k1xθ1

1(k)

−l1c0(k)+ε12x2(k)

−x1(k)h+

i=1

α1ix1(k)√hηi

i=1

β1ix1+θ1

1(k)√hηi

i=1

γ1ix1(k)c0(k)√hηi,

x2(k+ 1) = x2(k) + x2(k)r2−k2xθ2

2(k)

−l2c0(k)+ε21x1(k)

−x2(k)h+

i=1

α2ix2(k)√hηi

i=1

α2ix1+θ2

2(k)√hηi

i=1

γ2ix2(k)c0(k)√hηi,

where ηirepresents independent Gaussian random

variables with a standard normal distribution, denoted

as N(0,1). In the following figures, we choose

c0(t) = 0.1+0.05 sin(t).

5.1 Extinction

Example 5.1 We choose x1(0) = 0.7,r1= 0.06 ,

l1= 1,k1= 0.7,α1= 0.5,β1= 0.95 ,θ1= 0.5,

ε12 = 0.9,γ1= 0.05,x2(0) = 0.3,r2= 0.05,

l2= 1,k2= 0.8,α2= 0.51 ,β2= 0.85,θ2= 0.6,

ε21 = 0.8and γ2= 0.1. This gives

M−1

2m2=−0.09125 <0.

Therefore, the extinction condition stated in Theorem

3.2 is fulfilled, as confirmed by the computer simula-

tions depicted in Figure 1.

0 500 1000 1500 2000 2500 3000 3500 4000

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

Path of x1(t) and x2(t)

Time(t)

x1(t)

x2(t)

Fig. 1. Trajectories of x1and x2for SDE (2).

5.2 Persistence and Stationary distribution

Example 5.2 Set x1(0) = 0.7,r1= 0.7,l1= 0.2,

k1= 0.4,α1= 0.15 ,γ1= 0.1,β1= 0 ,θ1=

1,ε12 = 0.35,x2(0) = 0.8,r2= 0.7,l2= 0.1,

k2= 0.3,α2= 0.2γ2= 0.1,β2= 0,θ2= 1 and

ε21 = 0.4. This gives

r1−l1−ε12 −1

2∥α1+γ1∥2≈0.2969 >0,

and

r2−l2−ε21 −1

2∥α2+γ2∥2≈0.5167 >0.

Hence, the persistence condition of Corollary 3.1 is

verified. The computer simulations in Figure 2and

Figure 3support this outcome.

WSEAS TRANSACTIONS on MATHEMATICS

DOI: 10.37394/23206.2023.22.67

A. Nait Brahim, B. Harchaoui, S. Boutouil,

M. El Idrissi, S. Aznague, A. Settati,

A. Lahrouz, M. El Jarroudi

E-ISSN: 2224-2880

617

Volume 22, 2023

0 500 1000 1500 2000 2500 3000 3500 4000

0.5

1.5

2.5

Path of x1(t)

Time(t)

x1(t)

tZt

x1(s)ds

Fig. 2. Trajectories of x1and x2for SDE (2).

0 500 1000 1500 2000 2500 3000 3500 4000

0.5

1.5

2.5

3.5

Path of x2(t)

Time(t)

x2(t)

tZt

x1(s)ds

Fig. 3. Trajectories of x1and x2for SDE (2).

Example 5.3 Set x1(0) = 0.7,r1= 0.4,l1= 0.1,

k1= 0.4,α1= 0.15,γ1= 0.05,β1= 0.15,θ1=

0.85 ,ε12 = 0.1,x2(0) = 0.8,r2= 0.5,l2= 0.15,

k2= 0.3,α2= 0.2,γ2= 0.1,β2= 0.3,θ2= 0.95

and ε21 = 0.15. This gives

r1−l1−ε12 −1

2∥α1+γ1∥2= 0.18 >0,

and

r2−l2−ε21 −1

2∥α2+γ2∥2= 0.155 >0.

Therefore, the condition for a stationary distribution

as stipulated in Theorem 4.2 is validated. This

assertion is substantiated by the computational

simulations depicted in Figure 4.

Fig. 4. Kernel density functions of (x1, x2).

6 Discussion

Investigating stochastic dynamics using the Gilpin-

Ayala model, particularly within dispersed polluted

environments, is essential in modern ecological re-

search. This investigation delves into the intri-

cate interplay between environmental processes and

stochastic fluctuations in scenarios where pollutants

are disseminated throughout ecosystems. The Gilpin-

Ayala model is a fundamental tool for understanding

population dynamics, offering insights into species

interactions and coexistence dynamics. However, ex-

tending this model to encompass stochastic elements

in polluted environments introduces a more realistic

depiction of ecological systems, where inherent ran-

domness and external perturbations play pivotal roles.

One of the primary implications of incorporating

stochasticity in the Gilpin-Ayala model within pol-

luted environments is its capacity to capture the vari-

ability and uncertainty inherent in real-world ecolog-

ical systems. As an exogenous factor, pollution intro-

duces fluctuations in vital parameters such as growth

rates, mortality rates, and species interactions. Con-

sequently, the model’s stochastic version enables ex-

ploring how these uncertainties influence species sur-

vival, coexistence, and potential extinction. By study-

ing the stochastic dynamics of the Gilpin-Ayala model

in dispersed polluted environments, researchers can

attain more profound insights into the resilience of

ecosystems against pollution-induced disturbances.

Examining critical factors influencing extinction and

persistence in such scenarios provides a comprehen-

sive understanding of how ecosystems respond to en-

vironmental challenges. Moreover, this investiga-

tion advances theoretical ecology and applied envi-

ronmental science. The development of mathemati-

cal frameworks to analyze stochastic ecological mod-

els underlines the interdisciplinary nature of this re-

search, fostering collaboration between ecologists,

mathematicians, and statisticians. These models can

WSEAS TRANSACTIONS on MATHEMATICS

DOI: 10.37394/23206.2023.22.67

A. Nait Brahim, B. Harchaoui, S. Boutouil,

M. El Idrissi, S. Aznague, A. Settati,

A. Lahrouz, M. El Jarroudi

E-ISSN: 2224-2880

618

Volume 22, 2023

aid in devising strategies for mitigating the adverse ef-

fects of pollution on ecosystems, informing conserva-

tion efforts, and facilitating sustainable resource man-

agement. In conclusion, delving into the stochastic

dynamics of the Gilpin-Ayala model in dispersed, pol-

luted environments opens up new avenues for under-

standing the intricate dynamics of ecological systems

in the face of uncertainty and external stressors. This

endeavor enhances our theoretical understanding of

ecology and offers practical implications for manag-

ing and conserving biodiversity in polluted environ-

ments.

7 Conclusion

In conclusion, investigating the stochastic dynamics

of the Gilpin-Ayala model in dispersed polluted en-

vironments has provided valuable insights into the

complex behavior of species populations and pollu-

tant concentrations. The incorporation of stochastic

elements has allowed for a more realistic represen-

tation of the inherent variability and unpredictability

in these systems. The findings have underscored the

importance of considering stochastic dynamics when

studying and managing dispersed polluted environ-

ments, contributing to a more comprehensive under-

standing of their ecological dynamics and develop-

ing effective environmental conservation and pollu-

tion control strategies.

References:

[1] R. M. May, Stability and Complexity in Model Ecosystems,

USA, Princeton University Press (2001).

[2] M. Liu, K. Wang, Persistence and extinction of a stochastic

single species model under regime switching in a polluted

environment, Journal of Theoretical Biology, vol. 264, no.

3, pp. 934–944 (2010).

[3] M. E. Gilpin, F. G. Ayala, Global models of growth and

competition, Proceedings of the National Academy of Sci-

ences of the United States of America, vol. 70, no. 3, pp.

3590–3593 (1973).

[4] T. Caraballo, A. Settati, M. El Fatini, A. Lahrouz, A. Imlahi,

Global stability and positive recurrence of a stochastic SIS

model with Lévy noise perturbation, Physica A: Statistical

Mechanics and Its Applications, 523, 677-690 (2019).

[5] A. Settati, A. Lahrouz, A. Assadouq, M. El Fatini, M. El

Jarroudi, K. Wang, The impact of nonlinear relapse and

reinfection to derive a stochastic threshold for SIRI epi-

demic model, Chaos, Solitons & Fractals, 137, 109897,

https://doi.org/10.1016/j.chaos.2020.109897 (2020).

[6] M. El Idrissi, B. Harchaoui, A. N. Brahim, I. Bouzalmat, A.

Settati, A. Lahrouz, A sufficient condition for extinction and

stability of a stochastic SIS model with random perturba-

tion, WSEAS Transactions on Systems, 21, 367–371, DOI:

10.37394/23202.2022.21.40 (2022).

[7] A. Settati, A. Lahrouz, M. Zahri, A. Tridane, M. El Fatini,

H. El Mahjour, M. Seaid, A stochastic threshold to predict

extinction and persistence of an epidemic SIRS system with

a general incidence rate, Chaos, Solitons & Fractals, 144,

110690 (2021).

[8] S. Aznague, M. El Idrissi, A. N. Brahim, B. Harchaoui,

S. Boutouil, A. Settati, A. Lahrouz, M. El Merzguioui,

J. El Amrani, A Probabilistic SIRI Epidemic Model In-

corporating Incidence Capping and Logistic Population

Expansion, Appl. Math. Inf. Sci. 17, No. 5, 773–789,

http://dx.doi.org/10.18576/amis/170505 (2023).

[9] A. El Haitami, A. Settati, A. Lahrouz, M. El Idrissi, M. El

Merzguioui, A stochastic switched SIRI epidemic model in-

tegrating nonlinear relapse phenomena, Int. J. Dynam. Con-

trol, https://doi.org/10.1007/s40435-023-01256-9 (2023).

[10] T. C. Gard, Stochastic models for toxicant-stressed popula-

tions, Bulletin of Mathematical Biology, vol. 54, no. 5, pp.

827–837 (1992).

[11] M. Liu, K. Wang, Dynamics of a non-autonomous stochastic

Gilpin-Ayala model, Journal of Applied Mathematics and

Computing, vol. 43, no. 1-2, pp. 351-368 (2013).

[12] M. Liu, K. Wang, Survival analysis of stochastic single-

species population models in polluted environments, Eco-

logical Modelling, vol. 220, no. 9, pp. 1347–1357 (2009).

[13] Z. Geng, M. Liu, Analysis of Stochastic Gilpin-Ayala Model

in Polluted Environments, IAENG International Journal of

Applied Mathematics, 45 (2) (2015).

[14] G. Strang, Linear Algebra and its Applications, Thomson

Learning, Inc, London (1988).

[15] M. Liu, K. Wang, Survival analysis of a stochastic coopera-

tion system in a polluted environment, Journal of Biological

Systems, vol. 19, no. 2, pp. 183-204 (2011).

[16] M. Liu, K. Wang, Q. Wu, Survival analysis of stochastic

competitive models in a polluted environment and stochas-

tic competitive exclusion principle, Bulletin of Mathemati-

cal Biology, vol. 73, no. 9, pp. 1969–2012 (2011).

[17] J. Dhar, K. S. Jatav, Mathematical analysis of a delayed

stage-structured predator-prey model with impulsive diffu-

sion between two predator territories, Ecol. Complex. 16,

59–67 (2013).

[18] L. J. S. Allen, Persistence, Extinction, and Critical Patch

Number for Island Populations, J. Math. Biol. 24, 617–625

(1987).

[19] X. Zou, D. Fan, K. Wang, Effects of Dispersal for a Logistic

Growth Population in Random Environments, Abstr. Appl.

Anal, Article ID 912579 (2013).

[20] L. Arnold, Stochastic Differential Equations, Theory and

Applications, John Wiley and Sons, New York, NY, USA

(1972).

[21] A. Friedman, Stochastic Differential Equations and Appli-

cations, Vol. 2, Probability and Mathematical Statistics, vol.

28, Academic Press, New York, NY, USA (1976).

[22] A. LAHROUZ, A. SETTATI, A note on stochastic Gilpin-

Ayala population model with dispersal, Differential Equa-

tions and Dynamical Systems, vol. 25, no. 3, p. 417–430

(2017).

[23] C. Zhu, G. Yin, Asymptotic properties of hybrid diffusion

systems, SIAM J. Control Optim, 46, 1155–1179 (2007).

WSEAS TRANSACTIONS on MATHEMATICS

DOI: 10.37394/23206.2023.22.67

A. Nait Brahim, B. Harchaoui, S. Boutouil,

M. El Idrissi, S. Aznague, A. Settati,

A. Lahrouz, M. El Jarroudi

E-ISSN: 2224-2880

619

Volume 22, 2023

[24] M. Liu, K. Wang, Persistence and extinction of a single-

species population system in a polluted environment with

random perturbations and impulsive toxicant input, Chaos,

Solitons and Fractals, vol. 45, no. 12, pp. 1541–1550 (2012).

[25] X. Mao, Stochastic Differential Equations and Applications,

Horwood Publishing Limited, Chichester (1997).

[26] X. R. Mao, G. Marion, E. Renshaw, Environmental Brow-

nian noise suppresses explosions in populations dynamics,

Stochastic Processes and their Applications, vol. 97, no. 1,

pp. 95–110 (2002).

[27] E. P. Kloeden, Numerical Solution of Stochastic Differential

Equations, Springer, New York (1992).

Contribution of individual authors to

the creation of a scientific article

(ghostwriting policy)

A. Nait Brahim, B. Harchaoui, S. Boutouil: were

responsible for the conceptualization, validation, for-

mal analysis, writing - original draft, methodology,

writing - review & editing. M. El Idrissi, S. Az-

nague: have implemented the software, formal anal-

ysis, writing - original draft, writing - review & edit-

ing. A. Settati, A. Lahrouz, M. El Jarroudi: car-

ried out the validation, investigation, conceptualiza-

tion, writing - review & editing.

Sources of funding for research

presented in a scientific article or

scientific article itself

Firstly, we have been invited by your prestigious

mathematical journals to contribute to our research.

Furthermore, the authors affirm that they do not pos-

sess any discernible financial conflicts of interest

or personal affiliations that might influence the out-

comes and conclusions elucidated in this manuscript.

Creative Commons Attribution

License 4.0 (Attribution 4.0

International , CC BY 4.0)

This article is published under the terms of the Cre-

ative Commons Attribution License 4.0

WSEAS TRANSACTIONS on MATHEMATICS

DOI: 10.37394/23206.2023.22.67

A. Nait Brahim, B. Harchaoui, S. Boutouil,

M. El Idrissi, S. Aznague, A. Settati,

A. Lahrouz, M. El Jarroudi

E-ISSN: 2224-2880

620

Volume 22, 2023