Theory article

A stochastic Gilpin-Ayala mutualism model driven by mean-reverting OU process with Lévy jumps


  • Received: 28 November 2023 Revised: 26 January 2024 Accepted: 08 February 2024 Published: 23 February 2024
  • By using the Ornstein-Uhlenbeck (OU) process to simulate random disturbances in the environment, and considering the influence of jump noise, a stochastic Gilpin-Ayala mutualism model driven by mean-reverting OU process with Lévy jumps was established, and the asymptotic behaviors of the stochastic Gilpin-Ayala mutualism model were studied. First, the existence of the global solution of the stochastic Gilpin-Ayala mutualism model is proved by the appropriate Lyapunov function. Second, the moment boundedness of the solution of the stochastic Gilpin-Ayala mutualism model is discussed. Third, the existence of the stationary distribution of the solution of the stochastic Gilpin-Ayala mutualism model is obtained. Finally, the extinction of the stochastic Gilpin-Ayala mutualism model is proved. The theoretical results were verified by numerical simulations.

    Citation: Meng Gao, Xiaohui Ai. A stochastic Gilpin-Ayala mutualism model driven by mean-reverting OU process with Lévy jumps[J]. Mathematical Biosciences and Engineering, 2024, 21(3): 4117-4141. doi: 10.3934/mbe.2024182

    Related Papers:

    [1] Wenrui Li, Qimin Zhang, Meyer-Baese Anke, Ming Ye, Yan Li . Taylor approximation of the solution of age-dependent stochastic delay population equations with Ornstein-Uhlenbeck process and Poisson jumps. Mathematical Biosciences and Engineering, 2020, 17(3): 2650-2675. doi: 10.3934/mbe.2020145
    [2] Meili Li, Maoan Han, Chunhai Kou . The existence of positive periodic solutions of a generalized. Mathematical Biosciences and Engineering, 2008, 5(4): 803-812. doi: 10.3934/mbe.2008.5.803
    [3] Yongxin Gao, Shuyuan Yao . Persistence and extinction of a modified Leslie-Gower Holling-type Ⅱ predator-prey stochastic model in polluted environments with impulsive toxicant input. Mathematical Biosciences and Engineering, 2021, 18(4): 4894-4918. doi: 10.3934/mbe.2021249
    [4] Buyu Wen, Bing Liu, Qianqian Cui . Analysis of a stochastic SIB cholera model with saturation recovery rate and Ornstein-Uhlenbeck process. Mathematical Biosciences and Engineering, 2023, 20(7): 11644-11655. doi: 10.3934/mbe.2023517
    [5] Huili Wei, Wenhe Li . Dynamical behaviors of a Lotka-Volterra competition system with the Ornstein-Uhlenbeck process. Mathematical Biosciences and Engineering, 2023, 20(5): 7882-7904. doi: 10.3934/mbe.2023341
    [6] Xueqing He, Ming Liu, Xiaofeng Xu . Analysis of stochastic disease including predator-prey model with fear factor and Lévy jump. Mathematical Biosciences and Engineering, 2023, 20(2): 1750-1773. doi: 10.3934/mbe.2023080
    [7] Chun Lu, Bing Li, Limei Zhou, Liwei Zhang . Survival analysis of an impulsive stochastic delay logistic model with Lévy jumps. Mathematical Biosciences and Engineering, 2019, 16(5): 3251-3271. doi: 10.3934/mbe.2019162
    [8] Virginia Giorno, Serena Spina . On the return process with refractoriness for a non-homogeneous Ornstein-Uhlenbeck neuronal model. Mathematical Biosciences and Engineering, 2014, 11(2): 285-302. doi: 10.3934/mbe.2014.11.285
    [9] Giacomo Ascione, Enrica Pirozzi . On a stochastic neuronal model integrating correlated inputs. Mathematical Biosciences and Engineering, 2019, 16(5): 5206-5225. doi: 10.3934/mbe.2019260
    [10] Bingshuo Wang, Wei Li, Junfeng Zhao, Natasa Trisovic . Longtime evolution and stationary response of a stochastic tumor-immune system with resting T cells. Mathematical Biosciences and Engineering, 2024, 21(2): 2813-2834. doi: 10.3934/mbe.2024125
  • By using the Ornstein-Uhlenbeck (OU) process to simulate random disturbances in the environment, and considering the influence of jump noise, a stochastic Gilpin-Ayala mutualism model driven by mean-reverting OU process with Lévy jumps was established, and the asymptotic behaviors of the stochastic Gilpin-Ayala mutualism model were studied. First, the existence of the global solution of the stochastic Gilpin-Ayala mutualism model is proved by the appropriate Lyapunov function. Second, the moment boundedness of the solution of the stochastic Gilpin-Ayala mutualism model is discussed. Third, the existence of the stationary distribution of the solution of the stochastic Gilpin-Ayala mutualism model is obtained. Finally, the extinction of the stochastic Gilpin-Ayala mutualism model is proved. The theoretical results were verified by numerical simulations.



    As a common relationship among species, mutualism has been extensively studied by many experts and scholars. Mutualism models have also received a lot of attention in population dynamics [1,2,3]. For example, the Lotka-Volterra mutualism model, the most common model of interspecific relationships, has the following form [4]

    dxi(t)=xi(t)[riaiixi(t)+nj=1,jiaijxj(t)]dt,i=1,2,,n, (1.1)

    where xi(t) is the population size, ri is the intrinsic growth rate, aii>0 is the intraspecific competition coefficient, and aij>0(ji) is the effect of species j on species i. But, in the classical Lotka-Volterra mutualism model, the growth rate of each species is a linear function of the interacting species [5], which is unreasonable in real life. In order to describe the actual problem more accurately, Ayala and Gilpin et al. [5] proposed a nonlinear model in 1973

    dxi(t)=xi(t)[riaiixθii(t)+nj=1,jiaijxj(t)]dt,i=1,2,,n, (1.2)

    where θi denotes the positive parameter of the modified Lotka-Volterra mutualism model.

    However, in nature, no species is deterministic and will be affected by various environmental factors. To describe these random perturbations in the environment, we consider that the growth rate ri of species in model (1.2) is linearly disturbed by Gaussian white noise [6,7,8,9]

    ri(t)=ri+σidBi(t)dt,i=1,2,,n.

    For any time interval [0,t], let ˜ri(t) be the time average of ri(t). Then, we can get

    ˜ri(t):=1tt0ri(s)ds=ri+σiBi(t)tN(ri,σ2it),i=1,2,,n,

    where N(,) is the one-dimensional Gaussian distribution.

    However, it is unreasonable to use a linear function of Gaussian white noise to simulate random perturbations in real life [10]. Obviously, the variance of the average growth rate ˜ri tends to at t0+. This causes an unreasonable result that the stochastic fluctuations in the growth rate ri(t) can become very large in a small time interval [11]. Therefore, some scholars have begun to consider the use of mean-reverting Ornstein-Uhlenbeck process to simulate random perturbations, that is, the intrinsic growth rate ri of model (1.2) has the form [12,13]

    dri(t)=βi[ˉriri(t)]dt+σidBi(t),i=1,2,,n, (1.3)

    where βi is the reversion rate, σi is the intensity of environmental fluctuation, ˉri is the mean recovery level, and βi,σi>0. The mean reversion of ri(t) to the constant level ˉri when βi>0 can be inferred from (1.3): if ri(t) has diffused above ˉri at some time, then the coefficient of the dt drift term is negative, so ri(t) will tend to move downwards immediately after, with the reverse holding if ri(t) is below ˉri at some time [14,15].

    Further, we can get the solution of the OU process (1.3). First, by multiplying eβit on both sides of (1.3) and then sorting, we can get

    eβitdri(t)+βieβitri(t)dt=βiˉrieβitdt+σieβitdBi(t).

    Then,

    d(eβitri(t))=βiˉrieβitdt+σieβitdBi(t).

    Integrating from 0 to t on the both sides of above formula, we get

    eβitri(t)ri(0)=ˉri(eβit1)+t0σieβisdBi(s).

    Thus, we have

    ri(t)=ˉri+[ri(0)ˉri]eβit+σit0eβi(ts)dBi(s), (1.4)

    where ri(0) is the initial value of the Ornstein-Uhlenbeck process ri(t). Then, we can get the expectation and variance of ri(t) as follows:

    E[ri(t)]=ˉri+[ri(0)ˉri]eβit,Var[ri(t)]=σ2i2βi(1e2βit).

    Thus, ri(t) obeys the Gaussian distribution N(ˉri+[ri(0)ˉri]eβit,σ2i2βi(1e2βit)), and σit0eβi(ts)dBi(s) obeys the Gaussian distribution N(0,σ2i2βi(1e2βit)). From the mean of ri(t), it should be obvious to see the mean reversion feature: When ri(0) deviates from ˉri either upward or downward, the degree of deviation decays at the rate of eβit and approaches ˉri. When t+, the asymptotic mean and variance are ˉri and σ2i2βi, respectively, which can be understood as stationary, long-run equilibrium mean and variance.

    But, in real life, in addition to small environmental disturbances such as white noise, there are also sudden environmental disturbances that cause significant changes in the survival status of species [16], such as earthquakes, hurricanes, epidemics, and so on [17,18]. These phenomena cannot be described by white noise, and the introduction of Lévy jumps in the basic model is a reasonable way to describe these phenomena [17,18]. So, we construct the following stochastic Gilpin-Ayala mutualism model driven by the mean-reverting OU process with Lévy jumps,

    {dxi(t)=xi(t)[(ri(t)aiixθii(t)+nj=1,jiaijxj(t))dt+Zγi(z)N(dt,dz)]dri(t)=βi[ˉriri(t)]dt+σidBi(t),,i=1,2,,n, (1.5)

    where xi(t),i=1,2,,n is the left limit of xi(t), modified parameter θi1,i=1,2,,n, and Bi(t),i=1,2,,n are independent standard Brownian motions defined on the probability space (Ω,F,{F}t0,P). N is a Poisson counting measure with characteristic measure v with v(Z)<, and Z is a measurable subset of (0,). ˜N represents a compensating random measure of Poisson random measure N, defined as ˜N(dt,dz)=N(dt,dz)v(dz)dt. In order to satisfy the corresponding biological significance, we assume that for all zZ, the jump diffusion coefficients γi(z)>1,i=1,2,,n.

    The model studied in this paper is improved on the basis of the classical Lotka-Volterra model, which no longer assumes linear exponential growth of the population and uses the mean reversion OU process to simulate small perturbations in the environment. This is a more reasonable method than assuming that the population parameters are linearly disturbed by Gaussian white noise. Furthermore, we also take into account the sudden disturbance of the population, so we introduce Lévy jumps to construct the model (1.5) studied in this paper. As far as we know, there are relatively few studies on such models, so it is very meaningful to study the properties of model (1.5).

    For convenience, the following definitions are taken in this article:

    For the sequence cij(1i,jn), we let

    ˇc=max1i,jncij,ˆc=min1i,jncij.

    For a symmetric matrix A of order n, we define

    λ+max(A)=supxRn+,|x|=1xTAx.

    Assumption 2.1. For any k{1,2,...,n}, there exists a constant c>0, and the following inequalities hold:

    (1)Z[|ln(1+γk(z))|(ln(1+γk(z)))2]v(dz)<c,(2)Z|γk(z)|v(dz)<c,(3)Z|(1+γk(z))q1|v(dz)<c.

    Assumption 2.2. For matrix A=(0a12a1na210a2nan1an20), there is

    12λ+max(A+AT)<aii,i=1,2,,n.

    Remark 2.1. Assumption 2.1 indicates that the interference intensity of Lévy noise on the system should not be too large. Assumption 2.2 shows that although system (1.5) is a mutualism system, the intensity of intraspecific competition is still greater than the intensity of interactions between species. Otherwise, if the interference intensity of Lévy noise to the system is too large and the interaction intensity of species is greater than the intraspecific competition intensity, the solution of the system may explode in finite time.

    Theorem 2.1. If Assumptions 2.1 and 2.2 hold, for any initial value (x(0),r(0))=(x1(0),,xn(0),r1(0),,rn(0))Rn+×Rn, there exists a unique solution (x(t),r(t))=(x1(t),,xn(t),r1(t),,rn(t)) of model (1.5) on t0, and it remains in Rn+×Rn with probability one.

    Proof. Noting that all the coefficients of model (1.5) satisfy the local Lipschitz condition, for any initial value (x(0),r(0)), the system has a unique local solution (x(t),r(t)) on t[0,τe), where τe is the explosion time of the solution. Therefore, to prove the solution (x(t),r(t)) is global, it is needed to prove τe= with probability one only. Hence, we take a sufficiently large p0>0 such that each component of (x(0),er(0)) falls within [1p0,p0]. For each integer p0 greater than p, we define the stopping time

    τp=inf{t[0,τe):xi(t)(1p,p)oreri(t)(1p,p),forsomei=1,2,,n}. (2.1)

    Obviously, τp is monotonically increasing as p increases. For convenience, let τ=limpτp, then ττe holds with probability one. Therefore, if τ=, then τe=. In the following, we use proof by contradiction to prove τ=. Suppose τ= does not hold with probability one, then there exist constants T>0 and ε(0,1) such that P(τT)>ε. So, there exists p1p0 such that

    P(τpT)ε,forallpp1. (2.2)

    Defining a C2-function V on Rn+×Rn

    V(x(t),r(t))=ni=1(xi(t)1lnxi(t)+r4i(t)4).

    When xi>0, we have the inequality xi1lnxi,1in, so V is a nonnegative function.

    Using the Itˆo formula, we can get

    dV=LVdt+ni=1σir3idBi(t)+ni=1Z[xiγi(z)ln(1+γi(z))]˜N(dt,dz), (2.3)

    where

    LV=ni=1(xi1)(riaiixθii+nj=1,jiaijxj)+ni=1βir3i(ˉriri)+ni=132σ2ir2i+ni=1Z[xiγi(z)ln(1+γi(z))]v(dz). (2.4)

    Then, there exists a constant N>0 such that

    LVni=1aiixθi+1i+ni=1aiixθii+ni=1nj=1,jiaijxixj+ni=1rixini=1ri+ni=1βiˉrir3ini=1βir4i+ni=132σ2ir2i+ni=1Zxiγi(z)v(dz)ni=1Zln(1+γi(z))v(dz)ni=1aiixθi+1i+ni=1aiixθii+ni=112λ+max(A+AT)x2i+ni=1|ri|xi+ni=1|ri|+ni=1βiˉrir3ini=1βir4i+ni=132σ2ir2i+ni=1xiZ|γi(z)|v(dz)+ni=1Z|ln(1+γi(z))|v(dz)N. (2.5)

    Substituting Eq (2.5) into (2.3), we have

    dVNdt+ni=1σir3idBi(t)+ni=1Z[xiγi(z)ln(1+γi(z))]˜N(dt,dz). (2.6)

    Taking the integral from 0 to τpT on both sides of Eq (2.6) and taking the expectation, we obtain

    EV(x(τpT),r(τpT))V(x(0),r(0))+NE(τpT)V(x(0),r(0))+NT. (2.7)

    When pp1, let Ωp={τpT}. From Eq (2.2), we can obtain P(Ωp)ε, and from the definition of τp, for each ωΩp such that one of xi(τp,ω),eri(τp,ω)(i=1,2,,n) is equal to p or 1p so that V(x(τp,ω),r(τp,ω)) is not less than (p1lnp),(1p1+lnp), or 14(lnp)4, we have

    V(x(τp,ω),r(τp,ω))min{p1lnp,1p1+lnp,14(lnp)4}.

    According to Eq (2.7), we can get

    V(x(0),r(0))+NTE[IΩp(ω)V(x(τp,ω),r(τp,ω))]εmin{p1lnp,1p1+lnp,14(lnp)4},

    where IΩp(ω) represents the indicator function of Ωp. Let p. Then, >V((x(0),r(0))+NT=, and thus we have a contradiction. Therefore, τ= holds with probability one. Theorem 2.1 is proved.

    Assumption 3.1. For any q>0, there is

    nj=1,ji(qq+1aij+1q+1aji)aii<0,i=1,2,,n.

    Remark 3.1. Assumption 3.1 indicates that, in the mutualism system (1.5), for any species in the system, the intensity of intraspecific competition is greater than the sum of the weighted average of interspecific competition intensity, otherwise the system may not have a bounded qth moment.

    Theorem 3.1. If Assumptions 2.1 and 3.1 hold, for any initial value (x(0),r(0))=(x1(0),,xn(0),r1(0),,rn(0))Rn+×Rn, the solution (x(t),r(t))=(x1(t),,xn(t),r1(t),,rn(t)) of model (1.5) has the property that

    E[xi(t)]qκ(q),i=1,2,,n

    for any q>0, where κ(q) is a continuous function with respect to q. That is to say, the qth moment of the solution (x(t),r(t)) is bounded.

    Proof. For any q2, defining a nonnegative C2-function V : Rn+×RnR+

    V(x(t),r(t))=ni=1(xqi(t)q+r2qi(t)2q).

    Applying the Itˆo formula to the function V, we obtain

    dV=LVdt+ni=1σir2q1idBi(t)+ni=1Z((xi+xiγi(z))qqxqiq)˜N(dt,dz),

    where

    LV=ni=1xqi(riaiixθii+nj=1,jiaijxj)+ni=1βir2q1i(ˉriri)+ni=12q12σ2ir2q2i+ni=1Z((xi+xiγi(z))qqxqiq)v(dz).

    Then,

    LVni=1aiixθi+qi+ni=1nj=1,jiaij(qxq+1iq+1+xq+1jq+1)+ni=1|ri|xqi+ni=1βiˉrir2q1ini=1βir2qi+ni=12q12σ2ir2q2i+ni=1xqiqZ|(1+γi(z))q1|v(dz)=ni=1aiixθi+qi+ni=1nj=1,ji(qq+1aij+1q+1aji)xq+1i+ni=1|ri|xqi+ni=1βiˉrir2q1ini=1βir2qi+ni=12q12σ2ir2q2i+ni=1xqiqZ|(1+γi(z))q1|v(dz). (3.1)

    Let η=qmin{β1,β2,βn}. Using the Itˆo formula again, we have

    d(eηtV)=ηeηtVdt+eηtdV=ηeηtVdt+eηt(LVdt+ni=1Z((xi+xiγi(z))qqxqiq)˜N(dt,dz)+ni=1σir2q1idBi(t))=eηt(ηV+LV)dt+eηt(ni=1Z((xi+xiγi(z))qqxqiq)˜N(dt,dz)+ni=1σir2q1idBi(t)). (3.2)

    Integrating from 0 to t on both sides of Eq (3.2) and taking the expected value, we obtain

    E(eηtV)=V(x(0),r(0))+Et0eηs(ηV+LV)ds. (3.3)

    Combining this with Eq (3.1), we have

    ηV+LVni=1ηxqiq+ni=1ηr2qi2qni=1aiixθi+qi+ni=1nj=1,ji(qaijq+1+ajiq+1)xq+1i+ni=1|ri|xqi+ni=1βiˉrir2q1ini=1βir2qi+ni=12q12σ2ir2q2i+ni=1xqiqZ|(1+γi(z))q1|v(dz)sup(x,r)Rn+×Rn{ni=1ηxqiq+ni=1ηr2qi2qni=1aiixθi+qi+ni=1nj=1,ji(qaijq+1+ajiq+1)xq+1i+ni=1|ri|xqi+ni=1βiˉrir2q1ini=1βir2qi+ni=12q12σ2ir2q2i+ni=1xqiqZ|(1+γi(z))q1|v(dz)}:=κ1(q). (3.4)

    Substituting Eq (3.4) into (3.3), we get

    E(eηtV)V(x(0),r(0))+Et0eηsκ1(q)ds.

    Then,

    eηtEVV(x(0),r(0))+eηt1ηκ1(q).

    Further,

    lim suptE[xqi(t)]qlim suptEV(x(t),r(t))qlim supt(V(x(0),r(0))eηt+eηt1ηeηtκ1(q))=qκ1(q)η:=κ2(q),i=1,2,,n.

    This means E[xqi(t)]κ2(q),i=1,2,,n,t0,q2. According to Hölder's inequality, for any ˜q(0,2), we obtain

    E[x˜qi(t)](E[x2i(t)])˜q2(κ2(2))˜q2,i=1,2,,n.

    Let κ(q)=max{κ2(q),(κ2(2))˜q2}. Then,

    E[xqi(t)]κ(q),i=1,2,,n,q>0.

    Theorem 3.1 is proved.

    Remark 3.1. Similar to the proof of Theorem 3.1, we have E[ri(t)]2qQ(q),i=1,2,,n,q>0.

    In this section, we give sufficient conditions for the existence of the stationary distribution of the solution of model (1.5), which reflects the persistence of species over long periods of time and is an important asymptotic property of population development. Many scholars have also studied the stability of the system. For example, Shao [19,20] studied the asymptotic stability in the distribution of stochastic predator-prey system with S-type distributed time delays, regime switching, and Lévy jumps, and also studied the stationary distribution of predator-prey models with Beddington-DeAngelis function response and multiple delays in a stochastic environment, and used different methods to analyze the stability of the systems according to the different disturbances on the models; Liu et al. [21] gave sufficient conditions for the distribution stability of a two-prey one-predator model with Lévy jumps. Before giving the theorem of the existence of stationary distributions, we give several lemmas.

    Assumption 4.1. aiinj=1,jiaji>0,βi>1,i=1,2,,n.

    Remark 4.1. Assumption 4.1 shows that the impact of intraspecific competition intensity on population density is greater than the sum of the growing-promoting effects of other species on the species, and the reversion rate of the intrinsic growth rate under the interference of OU processes should not be too small. Otherwise, the system may not have a stationary distribution.

    Lemma 4.1. Let Xa(t)=(x1(t),,xn(t),r1(t),,rn(t)) and X˜a(t)=(˜x1(t),,˜xn(t),˜r1(t),,˜rn(t)) be solutions of model (1.5) with initial values of a=(x1(0),,xn(0),r1(0),,rn(0))D and ˜a=((˜x1(0),,˜xn(0),˜r1(0),,˜rn(0))D, where D is any compact subset of Rn+×Rn. If Assumptions 2.1 and 4.1 hold, then the following equation holds:

    limt+(E|x1(t)˜x1(t)|++E|xn(t)˜xn(t)|+E|r1(t)˜r1(t)|++E|rn(t)˜rn(t)|)=0,a.s..

    Proof. Defining a function W

    W=|lnx1ln˜x1|++|lnxnln˜xn|+|r1˜r1|++|rn˜rn|.

    Then, we obtain

    d+W=ni=1(sgn(xi˜xi)d(lnxiln˜xi)+sgn(ri˜ri)d(ri˜ri))=ni=1sgn(xi˜xi)[(ri˜ri)aii(xθii˜xθii)+nj=1,jiaij(xj˜xj)]dt+ni=1sgn(ri˜ri)[βi(ri˜ri)]dtni=1aii|xθii˜xθii|dt+ni=1nj=1,jiaij|xj˜xj|dtni=1(βi1)|ri˜ri|dt. (4.1)

    Taking the integral on both sides of Eq (4.1) and taking the expectation, we obtain

    EWW(0)ni=1aiit0E|xθii˜xθii|ds+ni=1nj=1,jiaijt0E|xj˜xj|dsni=1(βi1)t0E|ri˜ri|ds.

    Noting EW(t)0, we then have

    ni=1aiit0E|xθii˜xθii|dsni=1nj=1,jiaijt0E|xj˜xj|ds+ni=1(βi1)t0E|ri˜ri|dsW(0). (4.2)

    Let θi=1,i=1,2,,n. Then,

    ni=1(aiinj=1,jiaji)t0E|xi˜xi|ds+ni=1(βi1)t0E|ri˜ri|dsW(0).

    Thus, according Assumption 4.1, we have

    E|xi˜xi|L1[0,+),i=1,2,,n.

    Therefore, according (4.2), we get

    ni=1aiit0E|xθii˜xθii|ds+ni=1(βi1)t0E|ri˜ri|dsW(0)+ni=1nj=1,jiaijt0E|xj˜xj|ds+.

    Then, we have

    E|xθii˜xθii|L1[0,+),E|ri˜ri|L1[0,+),i=1,2,,n.

    According to model (1.5), there are

    E(xi(t))=x(0)+t0[E(ri(s)xi(s))E(aiixθi+1i(s))+nj=1,jiaijE(xi(s)xj(s))]ds+t0E(xi(s))Zγi(z)v(dz)ds,E(ri(t))=ri(0)+t0[E(βiˉri)E(βiri(s))]ds,i=1,2,,n.

    Therefore, E(xi(t)) and E(ri(t)),i=1,2,,n, are continuously differentiable. According to Theorem 3.1 and Remark 3.1, we have

    dE(xi(t))dt12E(x2i(t)+|ri(t)|2)+12nj=1,jiaijE(x2i(t)+x2j(t))+cE(xi(t))12(κ(2)+Q(1))+(n1)ˇaκ(2)+cκ(1),dE(ri(t))dtβi|ˉri|+βiE|ri(t)|βi|ˉri|+βiQ(1)12.

    So, E(xi(t)),E(ri(t)),i=1,2,,n, are uniformly continuous. According to the Barbalat lemma, it can be concluded that limt+E|xi˜xi|=0,limt+E|ri˜ri|=0,a.s., and therefore Lemma 4.1 is proven.

    Here, in order to prove the following lemma, we introduce the following symbols. Define B(Rn+×Rn) as the set of all probability measures on Rn+×Rn, and for any two measures p1,p2B, define the metric dH as

    dH(p1,p2)=suphH|Rn+×Rnh(x)p1(dx)Rn+×Rnh(x)p2(dx)|,

    where H={h:Rn+×RnR|h(x)h(y)||xy|,|h()|1}.

    Lemma 4.2. If Assumptions 2.1 and 4.1 hold, for any aRn+×Rn, {p(t,a,)t0} is the Cauchy sequence in the space B(Rn+×Rn) with metric dH.

    Proof. For any fixed aRn+×Rn, we only need to prove for any ε>0 that there is a T>0 such that

    dH(p(t+s,a,),p(t,a,))ε,tT,s>0.

    This is equivalent to prove

    suphH|Eh(Xa(t+s))Eh(Xa(t))|ε,tT,s>0. (4.3)

    For any hH,t,s>0, we have

    |Eh(Xa(t+s))Eh(Xa(t))|=|E[E(h(Xa(t+s))Fs]Eh(Xa(t))|=|Rn+×RnEh(Xz0(t))p(s,a,dz0)Eh(Xa(t))|Rn+×Rn|Eh(Xz0(t))Eh(Xa(t))|p(s,a,dz0)2p(s,a,ˉDcR)+ˉDR|Eh(Xz0(t))Eh(Xa(t))|×p(s,a,dz0), (4.4)

    where ˉDR={aRn+×Rn|a|R},ˉDcR=(Rn+×Rn)ˉDR. According to Chebyshev's inequality, the transition probability {p(t,a,dz0t0)} is compact, i.e., for any ε>0, there exists a compact subset D=D(ε,a) over Rn+×Rn such that p(t,a,D)1ε,t0, where R is sufficiently large and we have

    p(s,a,ˉDcR)<ε4,s0. (4.5)

    According to Lemma 4.1, there exists T>0 such that

    suphH|Eh(Xz0(t))Eh(Xa(t))|<ε2,t>T,z0ˉDR. (4.6)

    Substituting Eqs (4.5) and (4.6) into (4.4), we have

    |Eh(Xa(t+s))Eh(Xa(t))|<ε,tT,s>0. (4.7)

    Since h is arbitrary, inequality (4.3) holds.

    Lemma 4.3 [22]. Let M(t),t0, be a local martingale with initial value M(0)=0. If limt+ρM(t)<, then limt+M(t)t=0 where ρM(t)=t0dM,M(s)(1+s)2,t0, and M,M(t) is the quadratic variational process of M(t).

    Lemma 4.4. If Assumption 2.1 holds, the solutions of model (1.5) follow that

    lim suptlnxi(t)t0,i=1,2,,n,a.s.. (4.8)

    Proof. Defining a function W(t)=(ni=1xi(t))q=w(t)q,q1, using the Itˆo formula, we can get

    LW=q(ni=1xi(t))q1ni=1[rixiaiixθi+1i+nj=1,jiaijxixj]+ni=1xqiZ[(1+γi(z))q1]v(dz)qwq1(ni=1|ri|xi+ni=112λ+max(A+AT)x2i)+ni=1cxqiqwqni=1|ri|+qn12|λ+max(A+AT)|wq+1+ncwqni=1 q2q+1|ri|2q+1+n2q22q+1wq+12+qn12|λ+max(A+AT)|wq+1+ncwq.

    Let θ>0 be sufffciently small and satisfy mθt(m+1)θ,m=1,2,.... It follows that

    E[supmθt(m+1)θwq(t)]=E[wq(mθ)]+I,

    where

    I=E[supmθt(m+1)θ|tmθLWds|]E[supmθt(m+1)θ|tmθ(ni=1 q2q+1|ri|2q+1+2nq22q+1wq+12+qn2|λ+max(A+AT)|wq+1+ncwq)ds|]2nq22q+1E[(m+1)θmθwq+12(s)ds]+qn2|λ+max(A+AT)|E[(m+1)θmθwq+1(s)ds]+ncE[(m+1)θmθwq(s)ds]+ni=1q2q+1E[(m+1)θmθ|ri(s)|2q+1ds]2nq22q+1θE[supmθt(m+1)θwq+12(t)]+qnθ2|λ+max(A+AT)|E[supmθt(m+1)θwq+1(t)]+ncθE[supmθt(m+1)θwq(t)]+q2q+1θni=1E[supmθt(m+1)θ|ri(t)|2q+1].

    Choose θ sufffciently small such that I<h(q). Therefore,

    E[supmθt(m+1)θwq(t)]2h(q).

    Let ε be an arbitrary positive constant. Based on Chebyshev's inequality, it follows that

    P{supmθt(m+1)θwq(t)>(mθ)1+ε}2h(q)(mθ)1+ε,m=1,2,.

    By the Borel–Cantelli lemma, there exists an integer-valued random variable m0(ω) such that for almost all ωΩ, when mm0, we have

    supmθt(m+1)θwq(t)(mθ)1+ε.

    Hence, for almost all ωΩ, if mm0 and mθt(m+1)θ, we have

    lim suptlnwq(t)lntlim supt(1+ε)ln(mθ)ln(mθ).

    Letting ε0, we have

    lim suptlnwq(t)lnt1,a.s.,

    then,

    lim suptlnw(t)lnt1q,a.s..

    Thus,

    lim suptlnw(t)tlim suptlnw(t)lnt×lim suptlntt0,

    and it follows that

    lim suptlnxi(t)t0,i=1,2,,n,a.s..

    Lemma 4.5. If Assumption 2.1 holds, ˉri+lim supt1tt0Zln(1+γi(z))v(dz)ds>0,i=1,2,,n, then populations xi(t) are weak persistent, a.s..

    Proof. According to the definition of weak persistence, we need to prove lim suptxi(t)>0,i=1,2,,n. If the conclusion is not true, then P(U)>0, where U={ω:lim suptxi(t,ω)=0,i=1,2,,n}. Applying the Itˆo formula to lnxi(t) and integrating from 0 to t, we have

    lnxi(t)t=lnxi(0)t+1tt0(riaiixθii+nj=1,jiaijxj)ds+1tt0Zln(1+γi(z))v(dz)ds+Mi(t)t,(i=1,2,,n), (4.9)

    where

    Mi(t)=t0Zln(1+γi(z))˜N(ds,dz),i=1,2,,n.

    By Assumption 2.1,

    Mi,Mi(t)=t0Z[ln(1+γi(z))]2v(dz)ds<ct,i=1,2,,n.

    From Lemma 4.3, we obtain

    limtMi(t)t=0,i=1,2,,n.

    On the one hand, combining the strong law of large numbers [22] and the definition of the Ornstein–Uhlenbeck process, we have

    limt1tt0ri(s)ds=ˉri,i=1,2,,n.

    If for all ωU, lim suptxi(t,ω)=0,i=1,2,,n, combining with Eq (4.9) we have

    0lim suptlnxi(t,ω)t=ˉri+lim supt1tt0Zln(1+γi(z))v(dz)ds>0,i=1,2,,n.

    As this contradicts the assumption P(U)>0, then lim suptxi(t)>0,i=1,2,,n.

    Theorem 4.1. If Assumptions 2.1 and 4.1 hold, ˉri+lim supt1tt0Zln(1+γi(z))v(dz)ds>0,i=1,2,,n, and then model (1.5) has a unique ergodic stationary distribution.

    Proof. To prove Theorem 4.1, first prove that there is a probability measure η()B such that for any aRn+×Rn, the transition probability p(t,a,) for Xa(t) converges weakly to η().

    According to Proposition 2.5 [23], weak convergence of probability measures is the concept of a metric, i.e., p(t,a,) weakly converging to η() is equivalent to the existence of a metric d such that limt+d(p(t,a,),η())=0.

    So, we only need to prove that, for any aRn+×Rn, there is

    limt+dH(p(t,a,),η())=0.

    From Lemma 4.2, {p(t,0,t0)} is the Cauchy sequence in the space B(Rn+×Rn) of the metric dH. So, there is a unique η()B such that

    limt+dH(p(t,0,),η())=0.

    By Lemma 4.1 and the triangle inequality, we have

    limt+dH(p(t,a,),η())limt+[dH(p(t,a,),p(t,0,))+dH(p(t,0,),η())]=0.

    That is, the distribution of X(t) weakly converges to η.

    By the Kolmogorov-Chapman equation, we know that η is constant. From Corollary 3.4.3 [24], it follows that η is strongly mixed. From Theorem 3.2.6 [24], we know that η is ergodic.

    In this section, we give sufficient conditions for species extinction. For convenience, model (1.5) is written in matrix form as

    {dx(t)=diag(x1(t),x2(t),,xn(t))[(r(t)Sxθ(t)+Ax(t))dt+Zγ(z)N(dt,dz)]dr(t)=diag(β1,β2,,βn)[ˉrr(t)]dt+σdB(t), (5.1)

    where

    x(t)=(x1(t),x2(t),,xn(t))T,r(t)=(r1(t),r2(t),,rn(t))T,S=diag(a11,a22,,ann),xθ(t)=(xθ11(t),xθ22(t),,xθnn(t))T,A=(ajh)n×n(ajj=0),γ(z)=(γ1(t),γ2(t),,γn(t))T,ˉr=(ˉr1,ˉr2,,ˉrn)T,σ(z)=(σ1(t),σ2(t),,σn(t))T.

    Assumption 5.1. There exists a set of positive constants c1,c2,,cn such that

    λ+max(12(CA+ATC)CS)0

    holds, where C=diag(c1,c2,,cn).

    Remark 5.1. In Assumption 5.1, the introduction of the constant ci,i=1,2,,n, indicates that the intraspecific competition intensity of the i-th population and the interspecific interaction intensity of the i-th population to the other n1 species changes by ci times. If ci1, the intraspecific competition intensity and interspecific competition intensity increase by ci times; if ci<1, it is weakened by ci times. Assumption 5.1 means that, under the action of ci, the intraspecific competition intensity of each species is greater than the average of the action intensity of the species on other species and the action intensity of other species on the species. Otherwise, the population might not go extinct.

    Theorem 5.1. If Assumptions 2.1 and 5.1 hold, for any initial value (x(0),r(0))Rn+×Rn, the solution (x(t),r(t)) of system (5.1) has the property that

    lim suptln|x(t)|tmax1inˉri+max1in{aii(θi1)θθiθi1i}+lim supt1tt0Zln(1+ˇγ(z))v(dz)ds,a.s..

    In particular, if max1inˉri+max1in{aii(θi1)θθiθi1i}+lim supt1tt0Zln(1+ˇγ(z))v(dz)ds<0, it implies limt|x(t)|=0, and then x(t) is extinct, a.s..

    Proof. Define a Lyapunov function

    V(x)=cTx=ni=1cixi,xRn+,

    where c=(c1,c2,,cn)T.

    Applying the Itˆo formula, we can get

    dV(x)=xTC[r(t)Sxθ(t)+Ax(t)]dt+ZxTCγ(z)N(dt,dz).

    Using the Itˆo formula for lnV(x) again, we have

    dlnV(x)=1VxTC[r(t)Sxθ(t)+Ax(t)]dt+Z[ln(V(x)+xTCγ(z))lnV(x)]N(dt,dz)=1VxTC[r(t)Sxθ(t)+Sx(t)Sx(t)+Ax(t)]dt+Z[ln(V(x)+xTCγ(z))lnV(x)]N(dt,dz),

    where

    1VxTCr(t)max1inri(t),
    1VxTC[Sxθ(t)+Sx(t)]max1in{aii(θi1)θθiθi1i},

    where we use the fact aiixθii+aiixiaii(θi1)θθiθi1i,i=1,2,,n, and

    1VxTC[Sx(t)+Ax(t)]λ+max(12(CA+ATC)CS)|x(t)|ˆc0,
    Z[ln(V(x)+xTCγ(z))lnV(x)]N(dt,dz)Zln(1+ˇγ(z))N(dt,dz).

    Substituting the above four inequalities into dlnV(x), we get

    dlnV(x)max1inri(t)+max1in{aii(θi1)θθiθi1i}+Zln(1+ˇγ(z))v(dz)dt+Zln(1+ˇγ(z))˜N(dt,dz).

    Integrating from 0 to t, we have

    lnV(x(t))lnV(x(0))t0max1inri(s)ds+t0max1in{aii(θi1)θθiθi1i}ds+t0Zln(1+ˇγ(z))v(dz)ds+M(t), (5.2)

    where

    M(t)=t0Zln(1+ˇγ(z)˜N(ds,dz).

    By Assumption 2.1,

    M,M(t)=t0Z[ln(1+ˇγ(z))]2v(dz)ds<ct.

    From Lemma 4.3, we achieve

    limtM(t)t=0.

    On the one hand, combining the strong law of large numbers [22] and the definition of the Ornstein–Uhlenbeck process, we have

    limt1tt0ri(s)ds=ˉri,i=1,2,,n.

    Then,

    limt1tt0max1inri(s)dsmax1inlimt1tt0ri(s)ds=max1inˉri.

    According to Eq (5.2), we obtain

    lnV(x(t))lnV(x(0))t1tt0max1inri(s)ds+1tt0max1in{aii(θi1)θθiθi1i}ds+1tt0Zln(1+ˇγ(z))v(dz)ds+1tM(t). (5.3)

    Taking the upper limit on both sides of Eq (5.3), we get

    lim supt1tlnV(x(t))max1inˉri+max1in{aii(θi1)θθiθi1i}+lim supt1tt0Zln(1+ˇγ(z))v(dz)ds,a.s..

    When max1inˉri+max1in{aii(θi1)θθiθi1i}+lim supt1tt0Zln(1+ˇγ(z))v(dz)ds<0, it implies limt|x(t)|=0, then x(t) is extinct, a.s.. Theorem 5.1 is proved.

    Remark 5.1. Lemma 4.5 and Theorems 4.1 and 5.1 have very important biological explanations. From the theoretical results obtained, it can be seen that when ˉri+lim supt1tt0Zln(1+γi(z))v(dz)ds>0,i=1,2,,n, population xi(t),i=1,2,,n, will be weakly persistent, and if the parameters of model (1.5) satisfy the conditions of Assumption 4.1, the system has a stationary distribution, which indicates the persistence of population growth. When max1inˉri+max1in{aii(θi1)θθiθi1i}+lim supt1tt0Zln(1+ˇγ(z))v(dz)ds<0, and the parameters of model (1.5) satisfy the conditions of Assumption 5.1, population x(t)=(x1(t),,xn(t)) will be extinct. That is, for every 1in, when ˉri+lim supt1tt0Zln(1+γi(z))v(dz)ds<aii(θi1)θθiθi1i, population xi(t),i=1,2,,n, will be extinct. So, the survival and extinction of the biological population of model (1.5) completely depend on the value of ˉri+lim supt1tt0Zln(1+γi(z))v(dz)ds.

    Remark 5.2. In the following we analyze the effects of white noise simulated by the Ornstein-Uhlenbeck (OU) process on species survival and extinction. Since the OU process acts on the intrinsic growth rate ri,i=1,2,,n, if model (1.5) is not affected by jump noise, the model takes the following form:

    {dxi(t)=xi(t)[ri(t)aiixθii(t)+nj=1,jiaijxj(t)]dtdri(t)=βi[ˉriri(t)]dt+σidBi(t),,i=1,2,,n.

    Using a similar method as above, it can be proved that when ˉri>0,i=1,2,,n, populations xi(t),i=1,2,,n, are weakly persistent; when max1inˉri+max1in{aii(θi1)θθiθi1i}<0, population x(t)=(x1(t),,xn(t)) will be extinct. That is, when ˉri<aii(θi1)θθiθi1i,i=1,2,,n, populations xi(t),i=1,2,,n, are extinct. Thus, when the system is only disturbed by OU process, the survival and extinction of the population is only related to the value of the average growth rate ˉri,i=1,2,,n, of the population.

    When ˉri>0,i=1,2,,n, the species only disturbed by the OU process are weakly persistent. If the system is affected by jump noise and satisfies ˉri+lim supt1tt0Zln(1+γi(z))v(dz)ds<aii(θi1)θθiθi1i,i=1,2,,n, the species are extinct. When ˉri<aii(θi1)θθiθi1i,i=1,2,,n, the species that are only disturbed by the OU process are extinct, but if there are jump noises such that ˉri+lim supt1tt0Zln(1+γi(z))v(dz)ds>0,i=1,2,,n, the species are weakly persistent. Therefore, it can be obtained that jump noise can make the survival system extinct and the extinction system survive.

    Remark 5.3. In the following we analyze the effect of the jump diffusion coefficient γi(z),i=1,2,,n, on population survival and extinction. If γi(z)<0,i=1,2,,n, then lim supt1tt0Zln(1+γi(z))v(dz)ds<0,i=1,2,,n, means that jump noise could accelerate the extinction; if γi(z)>0,i=1,2,,n, then lim supt1tt0Zln(1+γi(z))v(dz)ds>0,i=1,2,,n, means that jump noise is beneficial to the survival of the population.

    In order to verify the above theoretical results on the stochastic Gilpin-Ayala mutualism model (1.5), we use the Euler-Maruyama method [25] and the R language, and select appropriate parameters for numerical verification. The combination of parameters is shown in Table 1, and the data is from [11,26,27,28,29]. Consider the following stochastic Gilpin-Ayala mutualism model for two populations:

    {dx1(t)=x1(t)[(r1(t)a11xθ11(t)+a12(t)x2(t))dt+Zγ1(z)N(dt,dz)]dx2(t)=x2(t)[(r2(t)a22xθ22(t)+a21(t)x1(t))dt+Zγ2(z)N(dt,dz)]dr1(t)=β1[ˉr1r1(t)]dt+σ1dB1(t)dr2(t)=β2[ˉr2r2(t)]dt+σ2dB2(t), (6.1)

    Example 6.1. Letting v(Z)=1, and take the initial value of model (6.1) as x1(0)=0.11,x2(0)=0.2,r1(0)=0.2,r2(0)=0.1, choosing the combination A1 as the parameter values of model (6.1), and using the R language for numerical simulation, Figure 1 is obtained. By calculating, we have

    12λ+max(A+AT)a110.150.5=0.35<0,12λ+max(A+AT)a220.150.4=0.25<0.

    Then, Assumption 2.2 is satisfied. According to Theorem 2.1, the global solution of the stochastic Gilpin-Ayala population model (6.1) exists.

    Table 1.  Several combinations of biological parameters of model (6.1).
    Combinations Value
    A1 a11=0.5,a12=0.2,a21=0.1,a22=0.4,θ1=1,θ2=1,γ1=0.4,γ2=0.2,β1=2,β2=2,ˉr1=0.3,ˉr2=0.2,σ1=0.5,σ2=0.3
    A2 a11=0.55,a12=0.22,a21=0.21,a22=0.46,θ1=1.2,θ2=1.5,γ1=0.4,γ2=0.2,β1=2,β2=2,ˉr1=0.3,ˉr2=0.2,σ1=0.5,σ2=0.3,q=2
    A3 a11=0.28,a12=0.12,a21=0.18,a22=0.26,θ1=1.3,θ2=2,γ1=0.25,γ2=0.2,β1=1.3,β2=1.3,ˉr1=0.3,ˉr2=0.3,σ1=0.6,σ2=0.7
    A4 a11=0.4,a12=0.16,a21=0.12,a22=0.5,θ1=2,θ2=2,γ1=0.1,γ2=0.2,β1=2,β2=2,ˉr1=0.35,ˉr2=0.3,σ1=0.5,σ2=0.3

     | Show Table
    DownLoad: CSV
    Figure 1.  Global solution of stochastic system (6.1) with stochastic noises (σ1,σ2)=(0.5,0.3): (a), (b) are the global solution of x1(t) and x2(t) in three cases; (c), (d) are the global solution of r1(t) and r2(t) in two cases. The relevant parameters are determined by the combination A1.

    The red lines in Figure 1(a), (b) represent the solutions of populations x1,x2 in a deterministic environment without any disturbance. It can be seen that the development trend of the population is a smooth curve, and the population will not explode due to the limitation of environmental resources. The blue lines in Figure 1(a), (b) show the variation trend of the populations x1,x2 whose growth rate is disturbed by the OU process. The green lines in Figure 1(a), (b) represent the global solution of the population under the disturbance of the OU process and Lévy noise, and since the jump noise values are both positive, it indicates that the jump noise plays a role in promoting the population growth. Combined with the figure, it can be found that, compared with the other two situations, the population number also increases significantly at the same time under the positive Lévy jump interference. Lévy jumps represent some disturbances in the environment that cause sudden changes in the survival condition of the population. For example, when t=16,t=22 in Figure 1(a), (b), we can also see that the population number changes suddenly, which indicates the effect of Lévy jumps on the population.

    The red lines in Figure 1(c), (d) represent intrinsic growth rates r1,r2, while the blue lines in Figure 1(c), (d) represent population growth rates disturbed by the OU process, indicating that the interference of random environmental factors will make the growth rate r1(t),r2(t) fluctuate randomly under the interference of the OU process.

    Example 6.2. Letting v(Z)=1, taking the initial value of model (6.1) as x1(0)=0.11,x2(0)=0.2,r1(0)=0.2,r2(0)=0.1, choosing the combination A2 as the parameter values of model (6.1), and using the R language for numerical simulation, Figure 2 is obtained. By calculating, we obtain

    (qq+1a12+1q+1a21)a110.34<0,(qq+1a21+1q+1a12)a220.25<0.

    Then, Assumption 3.1 is satisfied. The numerical simulation results show that E(xq1),E(xq2) are less than κ(q), so E(xq1)κ(q),E(xq2)κ(q),q>0 hold, and Theorem 3.1 is verified.

    Figure 2.  Moment boundedness of solution of stochastic system (6.1) with q=2. The relevant parameters are determined by the combination A2.

    From the biological point of view, since the environmental resources are limited, no biological population can grow indefinitely, so we hope that the system solution is ultimately bounded. In Figure 2, letting q=2, we have E(x21)κ(2),E(x22)κ(2), which indicates that the final second moment of the population is bounded, which conforms to the laws of survival in the real world.

    Example 6.3. Letting v(Z)=1, taking the initial value of model (6.1) as x1(0)=0.11,x2(0)=0.2,r1(0)=0.2,r2(0)=0.2, choosing the combination A3 as the parameter values of model (6.1), and using the R language for numerical simulation, Figure 3 is obtained. By calculating, we get

    a11a21=0.280.18=0.1>0,a22a12=0.260.12=0.14>0,ˉr1+lim supt1tt0Zln(1+γ1(z))v(dz)ds0.3+0.2230.523>0,ˉr2+lim supt1tt0Zln(1+γ2(z))v(dz)ds0.3+0.180.48>0.

    Then, Assumption 4.1 and the conditions of weak persistent are satisfied. Figure 3(a), (c) represent the solution of x1(t),x2(t), and Figure 3(b), (d) represent the histogram of the solution of x1(t),x2(t). According Theorem 4.1, model (6.1) has a stationary distribution η().

    Figure 3.  Existence of stationary distribution. Left-hand panels show the simulations of the solutions x1(t) and x2(t) of stochastic system (6.1). Right-hand panels show the frequency histograms of x1(t) and x2(t) of stochastic system (6.1).

    As can be seen from Figure 3(a), (c), the values of population x1(t) are mostly between 1.5–3, and the values of population x2(t) of are mostly between 1.3–2.5, mainly concentrated in the middle region. Figure 3(b), (d) is the frequency histogram of populations x1(t),x2(t), shows a trend that high in the middle and low at both ends, and obeys normal distribution approximately. This indicates that if Assumption 4.1 and ˉri+lim supt1tt0Zln(1+γi(z))v(dz)ds>0,i=1,2, hold, the populations will continue to grow steadily over time, the population size will not change dramatically, and the different populations of the system will coexist harmoniously.

    Example 6.4. Letting v(Z)=1, taking the initial value of model (6.1) as x1(0)=0.1,x2(0)=0.1,r1(0)=0.2,r2(0)=0.1, choosing the combination A4 as the parameter values of model (6.1), and using the R language for numerical simulation, Figure 4 is obtained.

    Figure 4.  Extinction of stochastic system (6.1) with γ1(z)=0.1,γ2(z)=0.2. Populations x1(t) and x2(t) are extinct in the two cases. The relevant parameters are determined by the combination A4.

    According to the selected parameters, matrix A is (00.160.120), matrix S is (0.4000.5), and taking C=IR2×2, then

    λ+max(12(CA+ATC)CS)0.30.

    Futher,

    max1inˉri+max1in{aii(θi1)θθiθi1i}+lim supt1tt0Zln(1+ˇγ(z))v(dz)ds0.007<0,

    and then Assumption 5.1 is satisfied. According to Theorem 5.1, the stochastic Gilpin-Ayala population model (6.1) is extinct.

    According Remark 5.2, when max1inˉri+max1in{aii(θi1)θθiθi1i}<0, the populations x1(t),x2(t) are extinct when the populations disturbed only by the OU process. The red lines in Figure 4(a), (b) show the populations x1(t),x2(t) whose growth rate is disturbed by the OU process. When populations are disturbed only by the OU process, populations x1(t),x2(t) are extinct at t=20. The green lines in Figure 4(a), (b) represent the global solution of the population under the disturbance of the OU process and Lévy noise, population x1(t) is extinct at t=30 and population x2(t) is extinct at t=45. In this example, we let γ1(z)=0.1,γ2(z)=0.2, and according Remark 5.3, this indicates that when the Lévy noise value is greater than 0, the population growth is promoted, and the positive Lévy noise will delay the extinction of the population.

    In this paper, we study the dynamic behaviors of a stochastic Gilpin-Ayala mutualism model (1.5) driven by the mean-reverting OU process with Lévy jumps. The existence and uniqueness of the global solution, the moment boundedness of the solution, the existence of the stationary distribution and extinction of the stochastic Gilpin-Ayala mutualism model (1.5) are proved and verified by numerical examples. The existence and uniqueness of the global solution and the moment boundedness of the solution show that, the population shows a fluctuating growth trend under the interference of various random factors, and for any q>0, populations xi(t) (i=1,2,,n) have bounded q-th moments. The existence of the stationary distribution and extinction of the solution show that when ˉri+lim supt1tt0Zln(1+γi(z))v(dz)ds>0,i=1,2,,n, model (1.5) has a stationary distribution η(), which indicates the persistence of population growth, and the populations x(t) will be extinct when the conditions given by the assumption are satisfied.

    However, in model (1.5), only the influence of the OU process and Lévy jumps on the survival of the population were considered. But, in the real world, there are many environmental factors that affect the population, such as rainfall, drought, seasonal changes, etc..These are the questions we will be working on in the future.

    The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

    This study was supported by National Natural Science Foundation of China (11401085), Heilongjiang Province Postdoctoral Funding Program (LBH-Q21059), Fundamental Research Projects of Chinese Central Universities (2572021DJ04).

    The authors declare there is no conflict of interest.



    [1] X. Abdurahman, Z. D. Teng, Persistence and extinction for general nonautonomous N-species Lotka-Volterra cooperative systems with delays, Stud. Appl. Math., 118 (2007), 17–43. https://doi.org/10.1111/j.1467-9590.2007.00362.x doi: 10.1111/j.1467-9590.2007.00362.x
    [2] B. S. Goh, Stability in models of mutualism, Am. Nat., 113 (1979), 261–275. https://doi.org/10.1086/283384 doi: 10.1086/283384
    [3] J. X. Pan, Z. Jin, Z. E. Ma, Thresholds of survival for an N-dimensional volterra mutual istic system in a polluted environment, J. Math. Anal. Appl., 252 (2000), 519–531. https://doi.org/10.1006/jmaa.2000.6853 doi: 10.1006/jmaa.2000.6853
    [4] L. Stone, The stability of mutualism, Nat. Commun., 11 (2020), 2648. https://doi.org/10.1038/s41467-020-16474-4 doi: 10.1038/s41467-020-16474-4
    [5] M. E. Gilpin, F. J. Ayala, Global models of growth and competition, Proc. Natl. Acad. Sci. U.S.A., 70 (1973), 3590–3593. https://doi.org/10.1073/pnas.70.12.3590 doi: 10.1073/pnas.70.12.3590
    [6] X. Y. Li, X. R. Mao, Population dynamical behavior of non-autonomous Lotka-Volterra competitive system with random perturbation, Discrete Contin. Dyn. Syst., 24 (2009), 523–545. https://doi.org/10.3934/dcds.2009.24.523 doi: 10.3934/dcds.2009.24.523
    [7] M. Liu, K. Wang, Population dynamical behavior of Lotka-Volterra cooperative systems with random perturbations, Discrete Contin. Dyn. Syst., 33 (2012), 2495–2522. https://doi.org/10.3934/dcds.2013.33.2495 doi: 10.3934/dcds.2013.33.2495
    [8] M. Liu, P. S. Mandal, Dynamical behavior of a one-prey two-predator model with random perturbations, Commun. Nonlinear Sci. Numer. Simul., 28 (2015), 123–137. https://doi.org/10.1016/j.cnsns.2015.04.010 doi: 10.1016/j.cnsns.2015.04.010
    [9] X. H. Zhang, K. Wang, Asymptotic behavior of stochastic Gilpin-Ayala mutualism model with jumps, Electron. J. Differ. Equations, 2013 (2013), 1–17.
    [10] W. R. Li, Q. M. Zhang, M. Anke, M. B. Ye, Y. Li, Taylor approximation of the solution of age-dependent stochastic delay population equations with Ornstein-Uhlenbeck process and Poisson jumps, Math. Biosci. Eng., 17 (2020), 2650-2675. https://doi.org/10.3934/mbe.2020145 doi: 10.3934/mbe.2020145
    [11] B. Q. Zhou, D. Q. Jiang, T. Hayat, Analysis of a stochastic population model with mean-reverting Ornstein-Uhlenbeck process and Allee effects, Commun. Nonlinear Sci. Numer. Simul., 111 (2022), 106450–106468. https://doi.org/10.1016/j.cnsns.2022.106450 doi: 10.1016/j.cnsns.2022.106450
    [12] E. Allen, Environmental variability and mean-reverting processes, Discrete Contin. Dyn. Syst. - Ser. B, 21 (2016), 2073–2089. https://doi.org/10.3934/dcdsb.2016037 doi: 10.3934/dcdsb.2016037
    [13] X. R. Mao, T. Aubrey, C. G. Yuan, Euler-Maruyama approximations in mean-reverting stochastic volatility model under regime switching, J. Appl. Math. Stochastic Anal., 2006 (2014), 5–12. https://doi.org/10.1155/JAMSA/2006/80967 doi: 10.1155/JAMSA/2006/80967
    [14] R. A. Maller, G. Müller, A. Szimayer, Ornstein–Uhlenbeck processes and extensions, Handb. Financ. Time Ser., 2009,421–437. https://doi.org/10.1007/978-3-540-71297-8_18 doi: 10.1007/978-3-540-71297-8_18
    [15] A. Melnikov, A Course of Stochastic Analysis, Springer: Switzerland, 2023. https://doi.org/10.1007/978-3-031-25326-3
    [16] J. H. Bao, X. R. Mao, G. Yin, C. G. Yuan, Competitive Lotka–Volterra population dynamics with jumps, Nonlinear Anal. Theory Methods Appl., 74 (2011), 6601–6616. https://doi.org/10.1016/j.na.2011.06.043 doi: 10.1016/j.na.2011.06.043
    [17] J. H. Bao, C. G. Yuan, Stochastic population dynamics driven by Lévy noise, J. Math. Anal. Appl., 391 (2012), 363–375. https://doi.org/10.1016/j.jmaa.2012.02.043 doi: 10.1016/j.jmaa.2012.02.043
    [18] X. H. Zhang, K. Wang, Stability analysis of a stochastic Gilpin–Ayala model driven by Lévy noise, Commun. Nonlinear Sci. Numer. Simul., 19 (2014), 1391–1399. https://doi.org/10.1016/j.cnsns.2013.09.013 doi: 10.1016/j.cnsns.2013.09.013
    [19] Y. F. Shao, W. L. Kong, A predator–prey model with beddington–DeAngelis functional response and multiple delays in deterministic and stochastic environments, Mathematics, 10 (2022), 3378. https://doi.org/10.3390/math10183378 doi: 10.3390/math10183378
    [20] Y. F. Shao, Dynamics and optimal harvesting of a stochastic predator-prey system with regime switching, S-type distributed time delays and Lévy jumps, AIMS Math., 7 (2021), 4068–4093. https://doi.org/10.3934/math.2022225 doi: 10.3934/math.2022225
    [21] M. Liu, C. Z. Bai, M. L. Deng, B. Du, Analysis of stochastic two-prey one-predator model with Lévy jumps, Physica A, 445 (2016), 176–188. https://doi.org/10.1016/j.physa.2015.10.066 doi: 10.1016/j.physa.2015.10.066
    [22] R. A. Lipster, A strong law of large numbers for local martingales, Stochastics, 3 (1980), 217–228. https://doi.org/10.1007/978-3-642-23280-0 doi: 10.1007/978-3-642-23280-0
    [23] N. Ikeda, S. Watanabe, Stochastic Differential Equations and Diffusion Processes, Japan: North-Holland Publishing Company, 1981.
    [24] G. D. Prato, J. Zabczyk, Ergodicity for Infinite Dimensional Systems, Cambridge: Cambridge University Press, 1996.
    [25] D. J. Higham, An algorithmic introduction to numerical simulation of stochastic differential equations, SIAM Rev., 43 (2001), 525–546. https://doi.org/10.1137/S0036144500378302 doi: 10.1137/S0036144500378302
    [26] X. F. Zhang, R. Yuan, A stochastic chemostat model with mean-reverting Ornstein–Uhlenbeck process and Monod-Haldane response function, Appl. Math. Comput., 394 (2021), 125833. https://doi.org/10.1016/j.amc.2020.125833 doi: 10.1016/j.amc.2020.125833
    [27] Y. L. Cai, J. J. Jiao, Z. J. Gui, Y. T. Liu, W. M. Wang, Environmental variability in a stochastic epidemic model, Appl. Math. Comput., 329 (2018), 210–226. https://doi.org/10.1016/j.amc.2018.02.009 doi: 10.1016/j.amc.2018.02.009
    [28] Y. Zhao, S. L. Yuan, J. L. Ma, Survival and stationary distribution analysis of a stochastic competitive model of three species in a polluted environment, Bull. Math. Biol., 77 (2015), 1285–1326. http://dx.doi.org/10.1007/s11538-015-0086-4 doi: 10.1007/s11538-015-0086-4
    [29] I. R. Geijzendorffer, W. van der Werf, F. J. J. A. Bianchi, R. P. O. Schulte, Sustained dynamic transience in a Lotka–Volterra competition model system for grassland species, Ecol. Modell., 222 (2011), 2817–2824. https://doi.org/10.1016/j.ecolmodel.2011.05.029 doi: 10.1016/j.ecolmodel.2011.05.029
  • This article has been cited by:

    1. Katarzyna Pichór, Ryszard Rudnicki, Stochastic models of population growth, 2025, 22, 1551-0018, 1, 10.3934/mbe.2025001
  • Reader Comments
  • © 2024 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(1115) PDF downloads(58) Cited by(1)

Figures and Tables

Figures(4)  /  Tables(1)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog